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FOREWORD 


This report constitutes interim documentation of efforts performed on 
Contracts NASl-15783 and NASl-15795. The purpose of this report is the 
presentation of the GIM code modifications and the latest flowfield calcula- 
tions performed with the General Interpolants Method (GIM) computer code 
for NASA- Langley Research Center. The work is comprised of code de- 
scription, algorithm development, turbulence modeling, reacting chemistry 
and interactive inputs in addition to flowfield calculations of various prob- 
lems of interest. Inquiries concerning this report should be directed to: 

Lawrence W. Spradley 

Lockheed-Huntsville Research & Engineering Center 

4800 Bradford Drive 

Huntsville, AL 35807 

Telephone (205) 837-1800, ext. 249 
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1. INTRODUCTION AND SUMMARY 


The General Interpolants Method (GIM) code was developed to analyze 
complex flow fields which defy solution by simple methods. The code uses 
numerical difference techniques to solve the full three-dimensional time- 
averaged Navier -Stokes equations in arbitrary geometric domains. The 
numerical analogs of the differential equations are derived by representing 
each flow variable with general interpolant functions. The point of departure 
then requires that a weighted integral of interpolants be zero over the flow 
domain. By choosing the weight functions to be the interpolants themselves, 
the GIM formulation can produce the classical implicit difference schemes. 
Choosing the weight functions to be orthogonal to the interpolant functions 
produces explicit finite difference type discrete analogs. By appropriate 
choice of constants in the weight functions, the GIM becomes analogous to 
standard finite difference schemes such as centered, backward, forward, 
windward and multi-step predictor -corrector schemes. The GIM analogs, 
however, are automatically produced for arbitrary geometric flow domains 
and hence is a general point of departure and provides flexibility in the choice 
of differencing schemes. 

The GIM computer code was originally written for the GDC 7600 machine. 
The first effort that was accomplished under Contract to NASA-Langley was 
the conversion and reprogramming of the code for the GDC -STAR (now termed 
CYBER 203) vector processor. The GIM-STAR code was then exercised for 
three-dimensional exhaust flows for application to Scramjet engine studies. 

The next sequential study in this computational fluid dynamics effort consisted 
of the development and application of a parabolized GIM algorithm, computa- 
tion of the flow, including spillage, in a model aircraft inlet and investigation 
of linearized block implicit schemes for GIM application. These tasks were 
accomplished under the two subject contracts through a cooperative effort of 
the Hypersonic Aerodynamics and Hypersonic Propulsion Branches. 



The most current effort, which is the subject of this report, is a con- 
tinuation of the GIM code development and application on the CYBER 203 
machine. Objectives of this effort include the following. 


• Complete the development of the quasi-parabolic scheme 
for solving the parabolized Navier -Stokes equations using 
spatial marching/relaxation. 

• Adapt and verify the quasi-parabolic algorithm for com- 
puting subsonic forebody flow fields. 

• Complete the adaptation of the quasi-parabolic code to 
hyperbolic flow field computation. 

• Continue the investigation and implementation of linearized 
block implicit schemes for the GIM elliptic and parabolic 
codes. 

• Complete incorporation of a two -equation turbulent kinetic 
energy transport model and a multiple length scale model. 

• Incorporate interactive computer technology for user- 
oriented input of the geometry description for the GIM 
code. 

• Complete incorporation of a hydrogen-air equilibrium 
chemistry model. 

• Incorporate a hydrogen-air and a hydrocar bon -air finite 
rate chemistry model and synthesize the model into a 
global reaction model with one or more equations. 

• Compute the flow field for several configurations including 
a missile fuselage, a wing/body case, a turbulent boundary 
layer in a model inlet and mixing/reaction of hydrogen and 
air in a duct. 


Certain of these stated objectives were accomplished in full and others 
partially. The plan of attack was to organize a set of overall tasks, some 
of which are interrelated, aimed at meeting the major objectives. This report 
is organized into sections with each section independently presenting details 
of the major tasks. The following is a summary of these sections: 

SECTION 2: THE GIM CODE - VERSION D 

Formulation, coding and check out of GIM parabolized Navier -Stokes 
code was completed. The major effort on this task was expended on obtaining 


reliable coding for the marching algorithms in vector FORTRAN language. 

A relatively large amount of logic was necessary to implement the elliptic, 
hyperbolic and parabolic code in one package. The GEOM module was ex- 
tensively modified to compute the geometry data one-plane-at-a-time for use 
in the spatial marchers. The INTEG integration module can now operate with 
a full elliptic flow field or in a spatial marching mode, either parabolic or 
hyperbolic. The GIMPLT plotting module was also upgraded to heuidle the 
output of the new code. The version which is now online at the Langley Center 
is termed Version D. The input data for each module of this version can be 
generated interactively on a remote terminal. Version D contains many 
features not reported in the previous user's guide. The following items 
are especially noteworthy, 

• Elliptic, parabolic, hyperbolic solvers operational 

• Interactive input and runstream generation available 

• Explicit finite difference algorithms available in all 
three solvers; MacCormack implicit (MI) scheme 
operational in elliptic solver 

• Algebraic eddy viscosity and two -equation TKE turbu- 
lence models coded in the elliptic and parabolic solvers. 

Algebraic (Baldwin- Lomax) model only is currently 
checked out. 

• Ideal gas option operational for one or two components. 

Equilibrium, complete reaction, model in code for hydrogen- 
air mixture. 

• Geometry module can generate grids for all three solvers. 

The GIM/Version D code is available on pernhanent files at NASA- 
Langley. 

SECTION 3: IMPLICIT ALGORITHM DEVELOPMENT 

This task was concentrated on developing and implementing a linearized 
block implicit algorithm for the GIM code. The specific algorithm chosen is 
termed the MacCormack implicit (Ref. 1.1). The two-dimensional version of the 
scheme presented in the MacCormack AIAA paper was coded using the GIM 
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methodology. Boundary conditions were not fully presented in the paper, thus 
development of general boundaries was required. Formulation and coding of 
a three-dimensional version of the MI scheme were also accomplished. Check 
out and investigation of the new algorithm were done on a simple Gouette flow 
problem. The following are among the major findings of this task. 

• The MI scheme is not readily vectorizable. 

• The overall numerical stability is governed by the 
streamwise CFL number = 1.0. This can still be 
many times the global CFL based on a tight normal- 
direction mesh. 

• The method appears to have good potential for com- 
puting high Reynolds number boundary layer problems. 

• With both the explicit and implicit algorithms avail- 
able in GIM, a large variety of problems can be attacked. 

SECTION 4: TURBULENCE MODEL IMPLEMENTATION 

This task consists of implementing a number of turbulence models in the 
GIM code and applying these models to a two-dimensional simulation of an 
aircraft inlet. The first type of model is termed "algebraic" and consists of 
the Baldwin-Lomax (B-L) eddy viscosity model. This scheme has been im- 
plemented, in general, in the GIM code elliptic and parabolic solvers. The 
user has at his option, no viscosity, artificial viscosity only, laminar constant 
or Sutherland's law values or turbulent (B-L) coefficients. Also coded in 
general format in GIM is a two-equation turbulent kinetic energy (TKE) model. 
Conservation law form differential equations are written for the turbulent 
kinetic energy and for the dissipation of TKE. The B-L algebraic model has 
been checked out for an incompressible flow over a flat plate at several 
Reynolds numbers. The TKE models have not been fully verified at the time 
of this writing. The same case solved with laminar and the B-L model is 
being repeated with the two equation TKE scheme. 
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SECTION 5: CHEMISTRY MODEL DEVELOPMENT 


To enable the computation of chemically reacting flows, the code was 
modified under this task to include the option of ideal gas or a chemical equi- 
librium reaction scheme. The basic premise of chemical equilibrium treat- 
ment in flow fields is that the local reaction rates are much faster than the 
flow residence time. The consequence is that the reactions proceed to com- 
pletion at each point within the flow field. It has been shown that for chemical 
equilibrium the thermochemistry computations can be uncoupled from the 
flowfield solution. The thermochemistry data can then be communicated to 
the flowfield code via one of two methods. First, the flow state variables 
(P,p, T), thermodynamics (y, M^), and species distribution can be generated 
a priori for an isentropic expansion process from a given stagnation condition 
and the results stored in tabular form for use via a table look-up procedure. 
Second, the chemical equilibrium properties can be computed as an integral 
part of the flowfield solution by using the equilibrium calculation (Ref. 1.2) as a 
subprogram. In either case the results are the same. A third choice is to 
specify a specific gas system, simplify the reaction model and code the equa- 
tions. Rather than solving a full set of species equations, a global specie 
conservation law is solved. In this task, we have accomplished the following 
in terms of chemistry models: 


• The simple "complete reaction" model (type 3 above) 
has been coded for the hydrogen-air system. The coding 
has been checked out on a simple case and found to be 
correct. This simple model is now being used to com- 
pute the mixing and reacting flow of air over a plate 
with hydrogen injection. 

• The full equilibrium model for any gas system has been 
previously developed for other codes. These subprograms 
are currently being added to GIM to allow equilibrium 
computations for any gas system. No results have been 
obtained at this writing. 

• The finite rate reaction model has been examined for in- 
clusion in GIM. The equations and solution scheme can 
be coded in general, with the specific reactions to be 
allowed and the rate data being an input to the code. This 
development has been previously completed and is shown 
in Section 5. No GIM/finite rate coding has yet been done. 
This will be the major emphasis during the coming year's 
program. 
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SECTION 6: INTERACTIVE INPUT MODULE 


A new module has been added to the GIM code. The purpose of the 
module is to render the input and operation of GIM much easier on the 
CYBER 203. This interactive program, entitled "RUNGIM," uses the CDC 
FORTRAN Extended (FTN) Version 4.6 language and CDC NOS. 1.3 CYBER 
control language (CCL). RUNGIM essentially replaces the report-format 
input guides in Ref. 1.3. Upon entering RUNGIM, via an interactive terminal, 
for example, questions are asked and the user must respond with answers 
and data. The interactive module performs the following jobs; 

• Supplies actual input data for the geometry and integration 
modules 

• Provides program updates to use as input for the CDC 
UPDATE processor 

• Supplies dynamic dimension data to set the size of GIM 
for a specific problem 

• Sets up runstr earns, including control cards, file sizes, 
etc., for executing a GIM run. 

The interactive module is currently operational at Langley for the GIM GEOM 
and INTEG modules Version D. Section 6 of this report gives details on exe- 
cuting RUNGIM and shows examples of its use. 

SECTION 7: FLOW OVER A WING-BODY CONFIGURATION 

In addition to the development aspects of this work, some effort was 
expended on applying the code to problems of interest at Langley. A number 
of cases were run in conduction with the development work and under a sepa- 
rate ’’applications" contract with NASA-Langley. One of the main pure calcu- 
lations associated with this cooperative effort is a "wing -body" problem. The 
inviscid, supersonic, three-dimensional flow was computed with the GIM hyper- 
bolic marching solver. The configuration consists of an ogive-cylinder fuse- 
lage with a wedge-shaped wing attached. The flow was assumed to be inviscid 
for this initial solution. The problem is currently being repeated with viscous 
effects. The wing-body interference flow field was computed with the GIM/QH 



code using approximately 20,000 grid points and required 210 seconds to 
solve the full three-dimensional field on the CYBER 203. The solution of the 
full flow field is shown in Section 7 as computer generated contour maps. The 
regions of wing-body interference are clearly seen. The solution for surface 
pressures are compared to measured data for this configuration. The agree- 
ment is reasonable in some regions and shows some deviations in others. 
The inviscid treatment of leading edges cind corners contributed to this 
disagreement. 

An important part of a research effort is to publish the findings and 
results of the studies. In addition to the NASA Contractor Reports, the results 
of the GIM code development and applications have been published at pro- 
fessional society meetings and in the open literature. During the past year, 
two papers were presented at AIAA meetings, both of which were published 
in the AIAA Journal , one paper in the NASA-Langley Grid Generation Con- 
ference and a presentation at the American Physical Society-Division of 
Fluid Dynamics Meeting. 

To further enhance the GIM code's capability and to utilize its current 
potential, the following items are recommended. 

Development 

• Formulate and code the MacCormack implicit method for the 
GIM forward marching solvers 

• Complete the full equilibrium chemistry model coding and 
checkout 

• Write a full finite rate model subprogram with specific 
reaction data as input 

• Synthesize a global hydrogen-air finite rate model 

• Incorporate triangular elements into the geometry package 
to facilitate easier modeling and to allow better grid con- 
trol 
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Applications 


• Compute parallel mixing of hydrogen and air in a duct, for 
equilibrium and finite rate 

• Determine the flow field., with spillage, in a three-dimensional 
inlet 

• Perform a viscous calculation with the GIM/QP solver over 
the wing/body configuration. 


An important part of the next year's effort should be directed toward 
assisting Langley personnel in running the GIM code. A course is planned 
to acquaint users at Langley with aspects of the code that have changed or 
are new since the previous user -orientation given last year. 


SECTION 1 REFERENCES 


1-1. MacCormack, R. W., "A Numerical Method for Solving the Equations 
of Compressible Viscous Flow," AIAA Paper 81-0110, January 1981. 

1-2. Svehla, R.A., and B. J. McBride, "FORTRAN IV Computer Program 
for Calculation of Thermodynamic and Transport Properties of Com- 
plex Chemical Systems," NASA TN D-7056, January 1973. 

1-3. Spradley, L.W., J. F. Stalnaker and A. W. Ratliff, "Hyperbolic/ 

Parabolic Development for the GIM/STAR Code," NASA CR-3369, 
Langley Research Center, Hampton, Va., December 1980. 
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2. THE GIM CODE - VERSION D 


The completion of the GIM code parabolized algorithms was a major 
part of this year's effort. Inclusion of elliptic, hyperbolic and parabolic 
solvers was accomplished. As a result of these updates to the code, a number 
of modifications were made to each module: 

• GEOM - Generates geometric boundaries of the problem 

and the finite difference grid 

• MATRIX - Performs matrix operations to obtain finite 

difference operators for the GEOM grid 

• INTEG - Integrates the flow field with user specified 

initial and boundary conditions 

• GIMPLT - Displays grid and/or flow field in contour maps 

• RUNGIM - Interactively produces a runstream and input 

data for the GIM code 

With the inclusion of the updated capability, the lastest version (as o| this 
writing) of the code is termed GIM/Version D. The following paragraphs 
summarize the coding changes made to achieve Version D. Charts showing 
the revised input formats are also given. Sections 3 through 7 of this report 
presents the technical aspects of each task in detail. 

2.1 PROGRAM MODIFICATIONS 

Several modifications were made to the CYBER 203 version of the GIM 
code geometry module (GEOM). Some of the modifications were necessary 
to remain compatible with Version D of the integration module, and some were 
made to increase user facility and reliability. 

The most significant modification to the geometry module involved 
changing the output file logic so that data could be output a line (2-D) or 
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plane (3-D) at a time as required by the QP version of the integrator module. 
This modification greatly reduced the storage requirements of the geometry 
module for QP-type analyses. In order to be able to compute in either the 
elliptic or QP modes the output logic for both the geometry data and the analog 
data had to be modified extensively. 

Another modification was to incorporate logic to count the number of 
boundary nodes in a given geometry. The number is used in dimensioning 
the integer module but is difficult to determine manually for complex geom- 
etries. 

The final set of modifications involved rework on some of the logic in 
the module that was too difficult for users to understand. Changes were made 
to the circular arc, conical arc and edge of revolution algorithms to increase 
their utility and reliability. 

The MATRIX module for the elliptic solver remains essentially intact 
from the previous version. The reading/writing of files has been changed 
for compatibility with GEOM D. At this writing, there are two matrix modules 
on file at the computer center. 

• MATRIX D - For use with elliptic solver 

• MATQP D - For use with spatial marching 

solvers 

The coding of the matrix operations for the two different solvers necessitated 
construction of separate modules. 

The major portion of the coding changes for GIM/D was made to the 
integration module, INTEG. Considerable logic is necessary to allow inte- 
gration in arbitrary three-dimensional regions. The INTEG D module 
contains: 
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• Elliptic, hyperbolic and parabolic solves 

• Two and three dimensions 



• Inviscid, laminar or turbulent flow 

• One or two -gas (ideal) or simplified equilibrium 
chemistry 

• Explicit finite difference or the new " MacCormack" 
implicit algorithm 

• Interactive input of data. 

The storage and dimensioning of INTEG was also reworked to make optimum 
use of the high speed memory of the CYBER 203. 

The GIMPLT graphics module required some extensive reworking to 
make it compatible with the spatial marching routines. First, the program 
was simplified by eliminating all unnecessary coding. The original plotting 
routine had several options which are currently unused and thus these were 
eliminated. Next, internal variable names corresponded to symbols which 
were also confusing. These variable names have been changed and are now 
representative of actual usage as it applies to the GIM code. Some external 
variable names used as input were changed for consistency and readability. 
However, the order of input and location of variables has remained the same 
as the STAR version of the plotter. This is consistent with published docu- 
mentation. Also, the subroutines have been resequenced to aid tracability 
of "GO TO" and "DO" statements. Then, comments were placed throughout 
the subroutines to aid users in locating regions of interest for debugging or 
modifying purposes. The program has been simplified, rewritten, and in- 
ternally documented for ease of usage. 

Plot documentation has been improved by the addition of subtitles for 
contour and velocity vector plots. The subtitles are located at the bottom 
of each plot frame to the right of the main title. 

The size of the model that GIMPLT can handle has increased dra- 
matically from about 5,000 nodes last year to over 28,000 nodes currently. 
The primary limitations are the core size of CDC 7600 and the unavailability 
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of the CY203 for plot runs. The alternative has been to rewrite the plot program 
or to do cumbersome external data file manipulation to reduce the number of 
nodes which can be plotted for a given run. Fortunately, it has been possible 
to develop a way for the plotter to handle models of up to 28,000 nodes through 
variable manipulation within the plotter. This approach maximizes the usage 
of available core. To further increase the plotter capacity would require 
having more core available since the plotter must have at least enough room 
to store the X, Y and Z positions for each node. This modification allows one 
version of the plotter to be used for both small and large models without ex- 
ternal data file manipulation for models up to 28,000 nodes. 

2.2 STATE OF THE CODE 

The GIM/D code is cataloged on permanent files at NASA-Langley. 

It is recommended at the current time that the interactive input module 
described in Section 6 be used for inputting the GIM code. For com- 
pleteness and for information on the required program inputs, the follow- 
ing charts are included here. These summary input guides can be used 
in conjunction with the previous GIM documentation to fully describe the 
codes operation: 

GIM Manual - NASA CR 3157 ; 1979 

Hyper bolic/Parabolic — NASA CR 3369 > 1980 

CHARTS: Input Guides 

2-1 DYNMATD 

2-2 DYNDIMD 

2-3 GEOMD 

2-4 MATRIX D 

2-5 MATQPD 

2-6 INTEGD 

2-7 Notes on INTEGD 

2-8 GIMPLTD 
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Cards/Format 


NX IDIM NSPEC NXMAX NMATE NSAV METHOD 

(715) 


NX 

IDIM 

NSPEC 

NXMAX 

NMATE 


NSAV 


total number of nodes (elliptic) 
number of nodes in one cross -plane (QP) 

dimensionality {= 2 or 3) 

number of special nodes 

0 for GEOMD runs 

maximum number of nodes in a cross -plane 

number of nodes to be checked for mating 
e.g., for a surface of 20 nodes mated to a surface of 50 nodes 
NMATE > 50. When in doubt set NMATE to the number of 
nodes in the largest zone. 

number of planes processed per record on the geometry files. 
The actual number of nodes processed per record is 


(NSAV + 1) *NXMAX 
2 *NXMAX 


METHOD = 
> 


0 for elliptic runs 
0 for QP runs 


elliptic 

QP 


Chart 2-1 - DYNMAT D Input Data 
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Cards/Format 


MN IDIM ISPEC NSP MNB METHOD IREC 

(715) 


MN 

IDIM 

ISPEC 

NSP 

MNB = 
METHOD = 
IREC 


total number of nodes (elliptic) 

number of nodes in one cross -plane (QP) 

dimensionality 

two-gas flag =0 1 ideal gas 

= 1 2 ideal gas 

number of special nodes 

(equal to NSPEC input to matrix module for QP) 
number of boundary terms (IB ^ 9) 

= 0 program sets MNB = MN 
< 2 elliptic 
> 2 parabolic 

number of records on the GEOM file 
= 1 for QP 


Chart 2-2 - DYNDIM D Input Data 



Card Type 


Parameter Llst/Format 


1 HEADER(I)?I = 1, 72 

(12A6) 

2 NZONES, IDIM, ISTEP, IMATRX, IMATE 

(515) 

3 IWRITE, LWRITE, NWRITE 

(315) 

4 KC(I), 1=1,6 

(6A5) 

5 NSECTS 

(15) 

6 MAPE(I), 1=1,12 

(1215) 

7 MAPS(I), 1=1,6 

(615) 

8 (IBWL(I), I = 1,6), ITRAIN 

(715) 

9 (NNOD(I), 1=1, 3), (ISTRCH(I), 1=1,3) 

(615) 

10 DIVPI(I), I = 1, 3 

(3E10.4) 

11 [AETA(J,I), I = 1, NNOD(J)] , J = 1, IDIM 

(8E10.4) 

12 [(AC(I,K, J), I = 1,8), K = 1,5], J = 1,4 or 12 

(8E10.4) 

13 [AS(I,J), I = 1,8], J = 1,6 

(8E10.4) 

14 (PT(I, J), I = 1,5), J = 1,4 or 12 

(8E10.4) 

15 f(PMAX(I, K, J), I = 1, 5), ETAMAX(K, J), K = 1,4], 
J = 1, 4 or 12 

(6E10.4) 

Chart 2-3 - GEOMD Input Guide 

See GIM documentation for explanation of FORTRAN symbols, 

NASA CR 3157 and CR 3369. 
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Card 

Parameters 

Format 

1 

NDX, NDY, NDZ, ISNOPT* 

(415) 

2 

KC(I), 1= 1,6 

(6A5) 

3 

N1,IC,NT 

(315) 


Chart 2-4 - MATRIX D Input Guide 


See GIM documentation for explanation of FORTRAN symbols, 
NASA CR 3157 and CR 3369. 
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Card Type 


Parameter List/Format 
* 

1 ITITLE(I), I = 1, 78 

(13A6) 

2 KCl(I), 1= 1,3, KC2(I) (1= 1,3) 

(6A5) 

3 NPL.T(I), MI2(I), MI3(I), NATT(I), L2(I), L3(I) 

(1=: l,NZONES+l) 

(615) 

4 INI, nC, INPP, IPRMAX 

(415) 


Cards 1 and 2 are as previously defined. 


Card 3; ZONE/PLANE Data 


NPLT(I) 

MIZ(I) 

MI3(I) 

NATT(I) 

L2(I),L3(I) 


Number of cross planes in Zone I 
Number of nodes in t )2 in Zone I 
Number of nodes in Tj ^ in Zone I 

Plane of attachment of Zone 1=0 (or blank) for 1=1 

Grid coordinates for first mated node in Zone 1=0 
(or blank) for 1=1 


EXAMPLE: 5x5 node cross plane mated to 7 x 7 node cross plane 
(x's are mated nodes) 



^2 


Chart 2-5 - MATQP D Input Guide 


See GIM documentation for explanation of FORTRAN symbols, 
NASA CR 3157 and CR 3369. 



This input is terminated with a -1 in Columns 4 and 5, 

Card 4; Analog Print Specifications 

- First node to be printed 
= 0 no nodes printed 
= -1 terminates input 

- Nodal increment 
= 0 only INlst node printed 

- Number of planes to which this specification applies 
IPRMAX - Maximum number of print nodes 

= 0 all nodes possible under this specification will be printed 

Notes; 1. All the INPP must sum up to the total number of planes in the problem 
2. If only a -1 appears in Columns 4 and 5 no analogs will be printed. 

EXAMPLE 

0 0 21 1 

2 2 5 10 10 

12 6 6 

-1 


INI 

nc 

INPP 


No analogs are printed for the first 21 planes; analogs for the 2nd, 4th, 
6th, . . . 20th nodes are printed for the next five planes; analogs for the 1st, 
3rd, 5th, . . . node are printed for the last six planes. 


Notes on MATQPD Output 

1. No output appears for plane 1 (the initial input plane) 

2. For plane 2 or plane NATT(I) +1 for which the number of nodes 
in a cross -plane in zone I is greater than that in Zone I-l (i.e., 
planes which require flowfield input), analogs are printed for 
nodes INI through 2*NX (where NX is the number of nodes in a 
cross -plane) incremented as specified. 

3. For all other planes the node number printed on the output is 
increased by NX since nodes 1 to NX are in the preceding plane. 
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Chart 2-5 - (Concluded) 



Card Type 
1 


Parameter List/Format 
ITITLE(I), I = 1, 80 (80 Al) 


2 


2a 

3 


3a 

3b 

*JL. 

4 '" 


5 


+ 


5a 


+ 


6 


7 

7a 

8 

9 


IDIM, METHOD, ITMAX, IPRNT, ITSAVE, ISTART, 
lOTYPE, lUNITS, ITSTRT, IVISC, IDIST, 

ISPEC, IDS, IBOUND, ITHERM 
(1515) 

IPSQ, INORM, ISBSM 
(315) 

INFOUT, IJUMPO, JJUMPO, NIOUT, NJOUT, ICALC, 
AMFLW 

(615, ElO.O) 

INFINL, IJUMPI, JJUMPI, NHN, NJIN, ICALC, OUTMFL 
(615, ElO.O) 

INFINL,, IJUMPI, JJUMPI, NUN, NJIN, INFOUT, IJUMPO, 
JJUMPI, NIOUT, NJOUT, ICALC 
(1115) 

lORD, IPCVG, SUMCHK, SLPCHK 
(215, 2E10.0) 

NPLT(I), MI2(I), MI3(I), NATT(I), L2(I), L3(I) 

I = l,NZONES + 1 
(615) 

NN, NNX, NDX, NNY, NDY, NNZ, NDZ, NPM, KZONES 
(915) 

KST, KNX, KDX, KNY, KDY, KNZ, KDZ 
(715) 

DTIME, DTFAC, INCDT 

(2E10.0, 15) 

REALMU, REALK, GAMSl, GAMS2, WMl, WM2, DK,RK 
(8E10.0) 

VF(I), 1=1,8 

(8E10.0) 

EMU, ELAM, ERHO, ESPEC 
(4E10.0) 

KC(I), 1= 1,6 

(6A5) 

NNPM(I), NCPM(I), (NNCPM(I, J), J= 1,5), ANGPM(I); 

1= 1,NPM 

(715, ElO.O) 


Chart 2-6 - INTEG D Input Guide 


See GIM documentation for explanation of FORTRAN symbols, 
NASA CR 3157 and CR 3369. 
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Card Type 


11+ 

NCT(I, J,K), PXPM(I, J,K), PYPM(I, 
K = 1,4; J = 1, NCPM(I); 1 = 1, NPM 
(15, 2E10.0) 

12 

RHOZ, PZ, ASTAR, NINC, A, B 
(3E10.0, 15, 2E10.0) 

13 

NJ, INC, NTOT, ITAN, ITYPE 
(515) 

14 

RI, UI, VI, WI, PI, CSI 
(6E10.0) 

15+ 

Nl, IC,NT 

(315) 

15a* 

INI, nc, INPP, IPRMAX 
(415) 

16"' 

NJ, INC, NTOT, ITAN, ITYPE 
(515) 

16 a"' 

RI, UI, VI, WI, PI, CSI 
(6E10.0) 

" QP Only 


+Elliptic Only 



Chart 2-6 - INTEG D Input Guide (Concluded) 


5 ^ 

See GIM documentation for explanation of FORTRAN symbols, 
NASA CR 3157 and CR 3369. 
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Card 2: 


IVISC 


ITHERM 


0 - Inviscid Ealer Equations 

1 - Numerical damping only 

2 - Laminar viscosity and damping 

3 - BaldwLn/Lomax turbulence model 
Thermal boundary conditions for no -slip walls 

0 - Free thermal boundaries 

1 - Adiabatic walls 

2 - Constant temperature walls 


Card 2a: 


IPSQ 

INORM 

ISBSM 


= SUMSQ print increment 

= N, SUMSQ’ s printed every Nth iteration 

= 0 SUMSQ* s normalized by 1st iteration values 

> 0 SUMSQ’ s not normalized 

= 0 SUMSQ’ s of conserved variables calculated 

> 0 SUMSQ’ s of primitive variables calculated 


Card 4: 


lORD 

IPCVG 

SUMCHK 

SLPCHK 


1 1st order backward 

2 2nd order backward 
Convergence check increment 

N, convergence checked at every Nth iteration 
SUMSQ tolerance 

Convergence assumed if all SUMSQ < tolerance 
SUMSQ slope tolerance 

Convergence assumed if all slopes of the SUMSQ 
curves are < SLPCHK 


Chart 2-7 - Notes on the New Input for INTEGD 
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Cards 4a; 


These are the same as for MATQPD input. This sequence 
must be terminated with a -1 in Cards 4 and 5. 

Card 6; 

INCDT <0 - Implicit technique 


Card 7a; (If IVISC < 1 omit Card 7a) 

Sutherland's law: 

VF(1) = Sutherland's Constant 
= 1.45804E-5 g/cm-sec 

= 7.30350E-7 lb /ft-sec (R 

VF{2) = T,. 

= 110.33 K 
= 198.6 R 

Baldwin/Lomax Turbulence Model 
VF(3) = CCP = 1.6 
VF{4) = CKLEB = 0.3 
VF(5) = (CWK)"^/^ = 2.0 
VF(6) = Pr.j, = 0.9 
VF{7) = CMUTM = 14.0 


Note on variable viscosity: 


The viscosity at any node is computed as 


^ ^T 

where 

= numerical damping coefficient 
= REALMU + VF(1) * + T^') 

Hj. = Turbulent viscosity. 
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Chart 2-7 - (Continued) 



The laminar thermal conductivity is calculated through the 
laminar Prandtl number — which is now entered as REALK. 


Card 15a ; QP print control 


INI 
lie 
INPP 
IPRMAX , 


as described for MATQPD 


Note on plane numbering; 

The plane number which appears in the QP integration is the 
number of integration planes. That is, one less than the actual 
geometric planes since the first plane is input. The total number 
of integrated planes, MAXPLN, is thus one less than the total 
number of planes. All of the INPP in the print specifications 
above must sum up to MAXPLN. 


Cards 16 & 16a; 


Input initial conditions for added zone input surfaces. These 
assume the same form as cards 13 and 14. Only card input 
and USERIP options are available for these. This input is 
read at the time that the plane is integrated. The information 
is stored at the end of the primitive variable vectors. As a 
result the first node to be initialized, NJ, is 

NJ = 2 * NN + 1 

where NN is the number of non-input nodes in the cross -plane. 


Chart 2-7 - (Concluded) 



Card Type 


Parameter List/Format 


1 ITITLE 

(A40) 

2 NN, ITRSTR, ITRBLK, KDIM, ISPC 

(515) 

3 GAMMA, FACTOR, RGAS, PO, TO, RHOO 

(6E10.0) 


Specs. 

S-1 NPLT, STITLE, IVIEW, ISYM, ITHETl, 

lAXISl, ITHET2, IAXIS2, IXTABL, lYTABL, 
VFAC 

(15, 5X, A20, 815, ElO.O) 

S-2 NTYPE, JO, IJUMP, JJUMP, NI, NJ, IPRNT, 

JOO 

(815) 


Grid 

G-1 'GRID,' lOPT, ICSCLE, NSPECS, (ISPEC(I), 

I-l, NSPECS) 

(A4, IX, 15, 25X, 215, 715) 


VVEC 

V-1 'VVEC,' lOPT, NITER, ICSCLE, NSPECS, 

(ISPEC(I), 1=1, NSPECS) 

(A4, IX, 215, 20X, 215, 715) 

(NITER = for QP) 


Chart 2-8 - GIMPLTD Input Guide 
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VVEC 


V-2 (ISPEC(I), 1=8, NSPECS) (if NSPECS>7) 

(45X, 715) 

I-l (ITER(I), 1=1, NITER) 

(1615) 


Contours 


C-1 


C-2 


I-l 


L-1 


ITYPE, IPOT, NITER, NC, ITABLE, INCR, 
ICSCLE, NSPECS, ISPEC(I), 1=1, NSPECS) 

(A4, IX, 515, 5X, 215, 715) 

(ISPEC(I), 1=1, NSPECS) (if NSPECS>7) 
(45X, 715) 

(ITER(I), 1=1, NITER) 

(1615) 

(CVAL(I), 1=1, NC) 

(8E10.0) 


(No I-l cards for QP) 


Description of Input Data 


Card 1 

ITITLE - This problem identification title appears on the bottom left 
of the plot frame for Grid, Contour, and Velocity Vector 
plots as one line of forty characters. 

Card 2 

NN - The total number of nodes in the model which is being plotted. 

NN must be consistent with the GEOMD and INTEGD runs 
which generated the grid and flow field. 


Chart 2-8 - (Continued) 
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ITRSTR - The iteration number of the first set of flow field informa- 
tion on the INTEGD File 22 which is being used as input to 
the GIMPLTD module. (= -1 for QP). 

ISPC - The two gas flag used in INTEGD 

ISPC = 0 for a single gas 

= 1 for two gases 


Card 3 

RGAS - The gas constant as input in INTEGD 
Card S-1 

STITLE - The plot specification title. This subtitle will appear on 
the bottom right of grid, contour and velocity vector 
plots as one line of twenty characters. 


Card S-2 - Element generation control parameter 


NTYPE 


JO 

IJUMP 

JJUMP 

NI 

NJ 

JOO 


0 a single element is input 

2 a network of two -node line connectors is generated 

3 a network of 4 -node elements is generated con- 
connecting two zones with different numbering 
schemes 

4 a network of 4 -node elements is generated 
Default NTYPE = 0 

For NTYPE = 3, JO is the first node in the first string of 
nodes . 

For NTYPE = 3, IJUMP is the nodal increment for the first 
string of nodes. 

For NTYPE = 3, JJUMP is the nodal increment for the 
second string of nodes. 

For NTYPE = 3, NI is the number of elements generated. 

For NTYPE = 3, NJ = 1 

For NTYPE = 3 only, JOO is the first node in the second 
string of nodes. 


Example 


For NTYPE = 3, a network of 4-node elements, are generated connecting 
two zones with different numbering schemes as illustrated on the follow- 
ing page. 
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Chart 2-8 - (Continued) 



Where Jl, J2, J3 and J4 are the nodes connected by the element that 
the algorithm generates for 1=1 through NI. 

J1 = JO + (I-l) * IJUMP 
J2 = JOO + (I-l) * JJUMP 
J3 = J2 + JJUMP 
J4 = J1 + IJUMP 

NTYPE = 3, JO = I, IJUMP = 1, JJUMP = 10, NI = 4, NJ = 1, JOO = 10 

1 2 3 4 5 

J = 1 

10 20 30 40 50 

1=1 1=2 1=3 1 = 4 

Note on calculation of core sizes 

The total core requirements for plotting including the program storage, 
system library, and all data arrays is given by the following formula. 

KFLjq = KMAXjq + 24000jq 


where. 


KMAX^q = NWRDS + 3 * NN + 1 

= 1000 + 3 * NN + 1 (NN = number of nodes) 

KMAX is the dimension of the array "A" in the GIMPLTD MAIN program 
which is normally set to 46001 and must be changed to plot a larger model. 
KMAX is also a program variable which must be reset in the IvlAIN program 
for a larger model. 

For KlvIAX = 46001, the maximum number of nodes in a model to be 
plotted is 

Ivfisj = KlViAX - NWRDS - 1 


Chart 2-8 - (Continued) 
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NN = - }P^9 . 7 . } = 15000 nodes 


The field length requirement for KMAX = 46,001 is 


KFLjq = 46 , 001 jq + 24,OOOjq 


= 20,001jq 

= 210,561g 


A field length of 330,000g (110,592 j^q) permits the following 
KlvlAXjQ = 110,592 ^q - 24,000jq = 86,592jq 
NN = 86,592 - lOQ . Q .. - , 1 . = 28,530 nodes. 


Chart 2-8 - (Concluded) 



3. IMPLICIT ALGORITHM DEVELOPMENT 


3.1 INTRODUCTION 

Numerical solution of the unsteady Navier -Stokes equations by explicit 
finite difference techniques has a number of disadvantages. The most serious 
one, from a practical engineering viewpoint, is the small time steps which 
are usually required to maintain stability. Computation of boundary layer 
flows at high Reynolds number, for example, requires a fine grid near solid 
boundaries, hence very small time steps and long computer run times. 

Numerical treatment of the steady state parabolic form of the Navier- 
Stokes equations face many of the same difficulties as the elliptic form. The 
spatial marching step size is constrained by the small grid spacing required 
to resolve boundary layers normal to a solid wall. Thus, marching downstream 
great distances can result in impractically long run times. 


One apparent solution for these difficulties is the use of implicit methods, 
some of which are unconditionally stable for any size time step/ spatial march- 
ing step. These schemes are not without problems of their own in terms of 
their practical use. Among the major difficulties are the following; 

1. Implicit finite differences, in general, lead to systems of 
nonlinear algebraic equations when applied to the Navier- 
Stokes equations. These must either be solved directly 
or linearized in some manner. 

Z. Multi -dimensional implicit methods lead to very large 
systems of simviltaneous algebraic equations. Even for 
linear systems, the efficient solution is not practical due 
to large size of the coefficient matrix. 

3. Fully implicit methods cannot be programmed as efficiently 
as explicit methods on advanced vectorized machines such 
as the CYBER 203. 
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In a recent paper, MacCormack (Ref. 3-1) presented a new implicit 
method which circumvents the first three of the above mentioned difficulties. 
The new method contains two stages. The first stage utilizes the explicit, 
second order, predictor -corrector finite difference method developed by 
MacCormack in 1969 (Ref. 3-2) which is widely used in many codes, including 
the GIM code. The second stage removes the explicit stability requirements 
by numerically transforming the finite difference equations into an implicit 
form. The resulting matrix equations are either upper or lower block bi- 
diagonal equations and are easily solved. Morever, the method preserves 
the conservation form of the Navier -Stokes equations so that shock capturing 
techniques can be employed. Because of these advantages along with its 
straightforward extension to three dimensions, the MacCormack implicit 
method was selected for implementation in the GIM code. 

3.2 TWO-DIMENSIONAL DEVELOPMENT 

The GIM code formulation of the implicit MacCormack method begins 
with the Navier -Stokes equations in conservation law form. In two dimensions 
and by neglecting body force terms and heat sources, these equations can be 
written as 


8U 9E 8F 

8t 8x 8y 


= 0 


U 


P 

pu 

pv 

P<? 
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E 


pu 

PUV - 

iP^+Dn - UT^ - - q^ 


F 


pv 

PVU - 

pV^ + P - Tyy 

(pg + P) V - - VTyy - qy 


where 



and where 


p = mass density ^ 

V = y-component of velocity p 

(? = total energy M 

X = bulk viscosity coefficient k 

T = temperature t 

X, y = space coordinates 


= x-component of velocity 
= pressure 

= kinematic viscosity coefficient 
= thermal conductivity 
= time coordinate 
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p = (y- 1 ) P 


(ideal gas law) 


The subsequent development of the GIM code formulation parallels that of 
MacCormack up to the presentation of the finite difference equations (Eq. (8) 
of Ref. 3-1). Here the finite difference equations are written in a form com- 
patible with the GIM code as follows: 


A E^^ A.f” 

Auf = - ( y + .4 - ) 

1 ' Ax Ay ' 


(3.1a) 


A. A I A , B 


) = AU? 

1 1 


(3.1b) 


= uf+At6u“'*'^ 

1 1 1 


(3.1c) 


AU^^+1 

1 


A A 

Ax ^ Ay 


(3.2a) 


(I + At 


A lAl 


• ) (I + At . ) 6U = AU- 


(3.2b) 


^ (6U.”'*'^ - SU.^'*'^) 

1 1 2 ' X X > 


(3.2c) 


n = time level 


= node number 


= (H^). 


explicit 


^n+1 _ 

X ~ \dt I. . 
' '1,1 


implicit 
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and where 


Af. 

ST » Ax * Ay Ay difference operators given by 


(for any property Q) 


^xQ- 
+ 1 


NCON 


Air = 5D ®l+j‘^i+j = Forward analog of (g) 

j = 0 ‘ 


A Q 


NCON 


~ -Q- ■ 

Ax I-J I-J 

j = 0 


— Backward analog of 


(«. 


(3.3) 


A Q. 
+ • 

“S? 


NCON 




j = o 


'i+j ~i+j 


^ Forward analog of 


»). 


A Q. 

= F bf .Q. . 
Ay 1-3 i-j 

j = 0 


= Backward analog of 


(If). 


The coefficients at , a.”, ht and h7 used above are the finite difference coef- 
ficients at node i obtained from the GIM interpolant scheme for the explicit 
MacCormack method on an arbitrary mesh with NCON connections to adjacent 
nodes. The Jacobian matrices] A | and |b| are given by 


1*1 = 


and B = 


S'^ S 

y By 


with S , and as defined in Ref. 3-1 and presented in Appendix 3. A 


at the end of this section. 


The solution of Eqs. (3.1) and (3.2) proceeds as outlined in Ref. 3-1 with 
two exceptions. First, the use of the general difference operators given by 
Eq. (3.3) produces a set of difference equations that is either upper or lower 
block triangular. These equations are solved directly via backward or forward 
substitution with little more computational work than for the block bidiagonal 
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system, which is a special case for rectangular meshes. The second devia- 
tion from the process described in Ref. 3-1 involves the implementation of 
boundary conditions and is discussed next. 

3.3 BOUNDARY CONDITION DEVELOPMENT 

The solution process in the published implicit MacCormack method 

involves forward and backward "sweeping" over the finite -difference mesh 

points up to but not including the boundary points. Boundary values are fixed 

in terms of U and/or AU. Note, however, that in the GIM code formulation, 

the explicit -implicit sequence, Eqs.(3.1a), (3.1b) and (3.2a) (3.2b), operate 

9TJ 

with and on instead of AU. Boundary values are set in terms 

of via constraints derived from a quasi -variational technique (Ref. 3-3). 
dt 

To accommodate these boundary conditions, the GIM code formulation employs 
a set of pseudo-mesh points around the periphery of the actual mesh, as 
illustrated below: 



• Real Mesh Points x Pseudo -Mesh Points 
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OTJ 

At these pseudo-mesh points, the value of is set to zero for all time; 

the value of U is not set nor indeed ever required. The GIM code explicit 

step, Eqs.(3.1a) and (3.2a), proceeds as usual, utilizing values at real 

mesh points only. The implicit step, Eqs.(3.1b) and (3.2b), now involves 

forward and backward sweeping on all real mesh points including the boundary 

9U 

points. Boundary conditions are imposed on after the explicit step, Eqs. 

(3.1a) and (3.2a), and again after the implicit step, Eqs. (3.1b) and (3.2b). The 

pseudo-nodes ate used in the implicit differencing scheme to ensure that the 

resulting difference matrix is a strictly upper or lower block triangular 

9U 

matrix. The boundary value of S 0 at the pseudo-nodes is used to drive 

9TJ 

the solution to the steady-state condition - 0 everywhere. 

The boundary condition treatment described above is easily extended to 
non -rectangular and non-orthogonal meshes in both two and three dimensions. 
This flexibility is inherent in the general difference operators, Eqs. (3.3), 
employed in the GIM code and permits consideration of a wide range of geo- 
metrical configurations in a very convenient manner. 

3.4 THREE-DIMENSIONAL DEVELOPMENT 

The three-dimensional implicit GIM code formulation is very similar 
to the two-dimensional formulation. The factored operator and Jacobian 
diagonalization approach used in the implicit MacCormack method allow a 
simple straightforward extension to three dimensions with only a relatively 
small increase in computational labor. The three-dimensional implicit GIM 
code equations are presented in Appendix 3.B at the end of this section. 

3.5 RESULTS AND DISCUSSIONS 

The implicit MacCormack method has been formulated and coded to be 
compatible with the GIM code. Both the two- and three-dimensional versions 
are available as an option in the current INTEG module. Only limited exper- 
ience with the method has been acquired thus far and this almost exclusively 
with simple two-dimensional problems on rectangular meshes. However, 
this has been sufficient to establish several important characteristics of the 
method. 
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3.5.1 Vec tori zation 


Like other implicit methods, the implicit MacCormack method does not 
appear to be amenable to vectorization for use on the STAR or other vector 
machines. Thus the GIM code implicit routines execute much more slowly 
than do the highly vectorized explicit routines. Some gain has been achieved 
by partially vectorizing the computation of the eigenvalue vectors and the 
components of the S matrices. However, the bulk of the implicit routines 
including the backward and forward mesh sweeps and imposition of boundary 
conditions remains in scalar coding. As a result, it appears that implicit 
runs will execute about three times slower per iteration than explicit 
runs. Some additional gain may yet be achieved through refinement of the 
implicit coding as more experience is acquired. 

3.5.2 Stability 

The published implicit MacCormack method of Ref. 3-1 is described as 

being stable for unbounded At. However, this has not been the case with the 

implicit GIM code formulation. Thus far, the upper bound on At for stability 

appears to be the value of At at which the flow field computation begins to 

become implicit in the streamwise or dominant flow direction, i*e., 

CFL , . =1.0. For a two-dimensional problem on a rectangular 

streamwise ^ ° 

mesh with the streamwise direction in the x-coordinate direction, this upper 
bound is given by Eq. (9) of Ref. 3-1 as 


At < 


1/2 

j 

( u + c)/Ax + (2j//pAx ) 


In no instance has an implicit GIM code calculation been able to exceed this 
upper -bound and remain stable. Even calculations which have been integrated 
in time to an apparent steady state quickly become unstable if At is increased 
beyond the CFL , . <1.0 limit. 
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While there is yet no explanation for the apparent descrepancy between 
the reported and observed stability of the naethod, several items are being 
investigated. Variations in the boundary condition implementation, modifica- 
tion of the eigenvalue terms, and permutations of the x-y operator sequence 
and forward-forward/backward-backward sequence are all being examined 
to determine their effect on the stability limit. 

An important point to note is that even if the ^^^streamwise — 
stability limit cannot be relaxed or circumvented, the implicit GIM code 
option will remain a very useful and valuable addition to the code. Many 
problems of current interest require a highly refined mesh in directions 
normal to the dominant flow direction in order to resolve boundary layers 
and other viscous phenomena. In these instances, the implicit GIM code 
option can be used to run with a time step many times larger than that im- 
posed by explicit stability requirements (CFL < 1.0) and yet remain within 
the At bound imposed by the ^^^g(.j.eamwise — limit. As an example, 
a simple two-dimensional Couette flow has been integrated using the implicit 
GIM code option with time steps as large as 18 times that required by explicit 
stability requirements and remained stable. It is envisioned that the gain 
over the explicit stability limit could be much larger for highly refined meshes. 

3.5.3 Recommendations 

While the vectorization possibilities and stability limit should continue 
to be investigated, the implicit GIM code option can be used to calculate 
several more difficult and realistic flow fields, including a shock wave- 
boundary layer interaction problem, a two-dimensional inlet flow with spill- 
age, a three-dimensional internal duct flow, and a turbulent boundary layer 
problem. In the near future, the implicit GIM code option can be coupled to 
the Quasi -Parabolic (QP) routines to allow efficient computation of the flow 
over wing-body and missile configurations. Other possibilities included 
coupling the implicit routines with the chemistry routines or the two-equation 
turbulent kinetic energy model. 
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Appendix 3. A 

TWO-DIMENSIONAL IMPLICIT GIM CODE EQUATIONS 


Governing Equations 


E 


F 


where 



9E 9F 
9x 9y 


= 0 


U 


P 

pu 

Pv 

PS 


pu 


pu + P - T 


XX 


puv - r 
^ xy 

(P^ + P) U 


UT - VT " <1 
' XX ' xy ^x 


pv 

pvu - 

PV^ + P - Tyy 

(pg + P) V - ar^y - VTyy ' qy 
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and where 


p = mass density 

V = y- component of velocity 

(? = total energy 

X = bulk viscosity coefficient 

T - temperature 

X, y = space coordinates 


= x-component of velocity 
= pressure 

= kinematic viscosity coefficient 
= thermal conductivity 
= time coordinate 


p = (r-i)p 


2^ 2 
u + V 


(ideal gas law) 


Finite Difference Equations 


a.f” 

= - ( — t , ^ + — t - ?■ . ) 

1 ' Ax Ay ' 


J ^ ■ ' , 'Y 

P ( (I-At 0(I-At— .)6U™ = 


At 

J* 1 1 


1 


A E. 
- 1 

Ax 


A F. 
- 1 


A A ^ |B| — TT 

(I+^‘ — • ) (1+ At — ) su“+l = Au."*^ 


^ (6uf« - eu.“+*) 
1 1 £ ' 1 1 ' 


40 



where 


n = time level 


i = node number 
n+1 


= {W-f. , , - (f )” . 

' /i, explicit ' '1,1 


implicit 


A A_ A A_ 

and where and are difference operators given by 

(for any property Q) 


. ^ NCON 

E 

Ax i-w i+j i+j 

j = 0 


^ Forward analog of 


(S). 


NCON 

E 

j = 0 


A Q. 

= E 


^ Backward analog of 


(a 


A,Q. 


NCON 


-+ 1 
Ay 

^ T-l 

II 

^i+j ^i+j 

A Q. 

NCON 


- 1 
Ay 

= E 

b. .Q. . 

i-j i-j 


3 = 0 



2 Forward analog of 


(f). 


^ Backward analog of 


+ - + 

with finite difference coefficients a. , a. , b, and b. • 

1 ' 1 ’ 1 1 


Jacobian Matrices 


= S‘" D. S 
X Ax 


B I = S“^ D„ S 
y B y 
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where 


S 

X 


S 

y 



0 0 -l/c^ 

pc 0 1 

0 10 
-pc 0 1 


0 0 -1/c^ 

1 0 0 

0 pc 1 

0 -Pc 1 


= VyP/p a = 1/2 (u 



10 0 

-ll/p 1/p 0 

-v/p 0 1/p 

ap -u3 -v|3 

1 0 0 

-u/p 1/p 0 

-v/p 0 1/p 

a|3 -u3 -vP 


^ + v^) 3 = y - 1 
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= 

max 

||u 

1 + 

2V 1 Ax . 

PAx ■ 2 At ' 

^A2 

= 

max 

{!“ 

+ c 


^A3 

= 

max 

{u 

I + 


A4 

= 

max 

{lu 

- cl 

1 + — 0 0 ' 
' ^ PAx 2 At ’ 

^B1 

= 

max 

{Iv 

1 + 


^B2 

= 

max -< 

[Ivl 

1 + 

^ '7 if’ 

^B3 

= 

max j 

IN 

+ C 1 

1 + ^ - r if- “-o} 

^B4 

= 

max 1 


" 
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Appendix 3.B 

THREE-DIMENSIONAL IMPLICIT GIM CODE EQUATIONS 
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and wher e 


P = mass density 
S = total energy 
P = pressure 

fl = kinematic viscosity coefficient 
k = thermal conductivity 


u = X- component of velocity 

w = z -component of velocity 

V = y-component of velocity 

X = bulk viscosity coefficient 

T = temperature 

coordinates 

(ideal gas law) 


x,y,z,t = space and time 


P = (y - 1) p 
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and where ^ ^ ^ ^ and ^ are difference operators 

given by (for any property Q) 
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4. TURBULENCE MODEL IMPLEMENTATION 


4.1 REVIEW 

The ability to compute actual time-dependent turbulent flow fields is 
beyond the range of current and foreseeable computer resources. It is im- 
possible to refine the computational mesh to the extent necessary to resolve 
the wide range of eddy length scales involved without exceeding computer 
speed and memory bounds. Such is also the case with any attempt to deter- 
mine the instantaneous behavior of the rapid turbulent fluctuations about the 
mean flowfield solution. As a result the effects of turbulence must be de- 
scribed by modeling. 

To avoid resolving every possible mode of turbulent flow, the first step 
in turbulence modeling is to time-average the Navier-Stokes equations. In 
this process the conserved variables are written as a sum of a time-average 
mean flow and a fluctuating part. Casting the Navier-Stokes equations in 
terms of these averaged variables results in a set of governing equations 
similar to the original set with some new terms referred to as the Reynolds 
stress and turbulent heat transfer terms. Rubesin and Rose (Ref. 4-1) use 
time -mass -averaged variables to arrive at a set of governing equations which 
are particularly useful in compressible flow calculations. 

The Reynolds stress and turbulent heat transfer terms (hereafter re- 
ferred to only as Reynolds stress terms) are now unknowns in the governing 
equations and closure of the system must be provided through empirical 
relationships or "turbulence models." Marvin (Ref. 4-2) provides an excel- 
lent review of turbulence modeling in compressible flows. There are two 
major classes of turbulence models available. The first of these models 
takes advantage of the Boussinesq hypothesis which states that the Reynolds 
stress terms in the Navier-Stokes equations can be modeled as the product 
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of the phenomenological **eddy” viscosity and the appropriately averaged local 
velocity gradient. These models then depend only on local mean flow quan« 
titles and empirically determined constants. The second class involves 
models that use fundamental transport equations to obtain more directly the 
spatial variation of the Reynolds* stress throughout the flow. These models 
relate the turbulent shear stress (turbulent eddy viscosity) to the turbulent 
kinetic energy and are called the TKE models. 

Eddy Viscosity Models; There are many algebraic eddy viscosity models 
available to directly relate the Reynolds stress terms to the mean flow condi- 
tions. One of the most widely used is that developed by Cebeci (Ref. 4-3). It 
is a two-layer model which employs a Prandtl-Van Driest mixing length formu- 
lation in the near-wall region and the Clauser wake formulation in the outer 
region (or in wakes). This model and, in general, most of the algebraic 
closures perform well for calculating attached turbulent boundary layers in 
zero or mild pressure gradients. 

The Cebeci-Smith model was modified by Shang and Hankey (Ref. 4-4) 
to include a relaxation length for the eddy viscosity. This model was found 
to give good predictions in regions of flow separation under the influences of 
strong adverse pressure gradients. 

Another eddy viscosity model in wide use is that of Baldwin and Lomax 
(Ref. 4-5) This model was developed for use in two- and three-dimensional 
Navier-Stokes codes and is applicable to separated flows in shear layers, 
wakes and wall boundary layers. This model is similar to that of Cebeci 
but it avoids the necessity of determining the outer edge of the boundary layer 
(or wake) by employing an approximation to the Clauser formulation in the 
outer region. This approximation amounts to the use of a vorticity distribu- 
tion function rather than the boundary layer thickness to determine the length 
scales. 

Turbulent Kinetic Energy (TKE) Models; The modeling of turbulent 
kinetic energy in compressible turbulent shear flows is discussed in detail 
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in Ref. 4-2. In these models the governing equations are supplemented by 
additional transport equations for turbulent kinetic energy and the various 
Reynolds stresses (Ref. 4-3). For computational efficiency, the most widely 
used models add only one or two extra transport equations. These are com- 
monly referred to as one-equation and two-equation models. 

One -Equation Models — In its original form, the one -equation model 
involves the solution of the turbulent kinetic energy equations. The values 
for the constants involved in the method were obtained through correlation 
of a variety of experimental data (Ref. 4-6). However, the model did not in- 
corporate the compressibility corrections described in Ref. 4-7. This cor- 
rection involves the modification of the dissipation constant. The introduction 
of the compressibility correlation provides a radical improvement in the 
ability of the one-equation model to predict flowfield development at super- 
sonic Mach numbers. This corrected correlation is described in Ref. 4-8. 
Rubesin (Ref. 4-9) has also devised a one-equation model for compressible 
flows from an extension of Glushko*s (Ref. 4-10) incompressible model. 

Two-Equation Models — There have been a number of two-equation tur- 
bulent kinetic energy models developed. The fundamental differences between 
the various two-equation models and the one-equation model just described 
are that in the two-equation models a transport equation is written for the 
turbulence energy dissipation length scale, (or equivalently, the turbulence 
energy dissipation rate) and the Prandtl-Kolmogorov relation is used to ob- 
tain the turbulent eddy viscosity. The most highly developed two-equation 
model is that given in Launder et al. (Ref. 4-11). The use of this so-called 
*'k-e'* model is advocated and described by Launder and Spalding (Ref. 4-12). 
However, this two-equation model does not include a compressibility cor- 
rection, and predictions made using this model for highly compressible flows 
are not accurate. Pan describes an extension of this model to three-dimensional 
steady compressible flows (Ref. 4-13). 

In general, TKE models provide more flexibility and give better pre- 
dictions of turbulent flows than do the algebraic models. However, adjustments 
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of the empirical constants in these models are required depending on flow 
conditions. The user must weigh the increased accuracy of the TKE models 
against the increase in required computation time and storage in order to 
choose between the methods. 

4.2 MODELS CURRENTLY IN GIM 

At the time of this writing the GIM code user may choose between the 
following turbulence models; 

• The Baldwin-Lomax algebraic eddy viscosity model for 
elliptic and QP calculations in the presence of one wall. 

• The two-equation TKE model for two- and three- 

dimensional elliptic and QP calculations. 

The Algebraic Model; The eddy viscosity model presently in the code 
is a variation of the Baldwin-Lomax model (Ref. 4-5). The turbulence effects 
are modeled by an eddy viscosity coefficient which enters into the flow 
equations through a modified total viscosity; 

/X = +/i^ , (4.1) 

where the laminar viscosity il can assume a constant value or be calculated 
from the temperature via Sutherland's law. The thermal conductivity is deter- 
mined through the input of a laminar and turbulent Prandtl number (Pr and 
Pr^, respectively) as 


K/Cp = ]U/Pr + jU^/Pr^ 

The eddy viscosity is given by a two-layer model; 




^^T^nner 

^/^T^outer 


(4.2) 


(4.3) 
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The description of the formulae used to calculate the inner and outer solu- 
tions is given in Ref. 4-5 and will not be repeated here. The only variation 
of the model in the code from that in Ref. 4-5 is in the definition of r (called 
y in the reference). In the GIM code r is the distance to the node point of 
interest in the computational grid from the node point on the solid wall which 
lies along the row of nodes including the point of interest. The symbol r^ 
denotes the least value of r for which the inner and outer viscosities are 
equal. Thus, the GIM code version and the model in Ref. 4-5 are equivalent 
if one set of grid lines is maintained perpendicular to the wall, i.e., orthog- 
onal grids. It is planned to extend the present model to multiple walls. 


The TKE Model; The differential equation model presently in the code 
is a time -dependent extension of that proposed by Pan (Ref. 4-13). This model 
requires the solution of two additional equations. The first one is for the 
turbulent kinetic energy 


k = -i (u*^ + v*^ + w*^) (4.4) 

2 

where u* , etc., are the mean square turbulent fluctuations of the velocity 
about the mean velocity profile. The other equation describes the transport 
of the rate of dissipation of turbulent kinetic energy, e. The variables k and 
e are related through an extension of the Prandtl-Kolmogorov relation for the 
turbulent eddy viscosity 


- C^P kVe 


(4.5) 


where C is an empirical constant. The governing equations for the flow field 
are given in Fig. 4-1. and empirical constants. This model is in 

the code and has been compiled but no test cases have been completed as of 
this writing. A few questions remain to be addressed on the use of this model. 
The exact treatment of the boundary conditions on k and e at solid walls is in 
question since the codes which use these models do not treat real solid walls 
but use extrapolated "wall functions." This boundary treatment has been known 
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Fig. 4-1 - TKE Model in the GIM Code 
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to yield poor results in the near -wall region where the model predictions of 
the turbulent length scale, , vary dramatically from experimental 

results (Ref. 4-14). The ability of the model to compute turbulence 
in the presence of pressure gradients strong enough to cause separation or 
in recirculating flow is in doubt. The greatest obstacle to proper boundary 
treatment at solid walls is a paucity of experimental information from which 
to determine the behavior of the dissipation at the wall. For compressible 
flows it is not known how much information about density fluctuations is 
obscured by mass-weighted averaging; nor, has the ability of models averaged 
in this way to compute flows with large density or pressure gradients been 
verified. 

The need for artificial damping to stabilize the time -iterative solution 
of these equations has also not been studied. The specific form of the damping 
terms needs to be derived. 

4.3 RESULTS OF COMPUTATION 

The two-dimensional spillage problem defined in Ref. 4.15 was used as 
a test case. 

The turbulent flow field over the leading edge of this 25-deg compression 
surface of the model inlet was calculated. Figure 4-2 shows the 840-node 
computational grid. The velocity vectors, pressure, and Mach contours are 
shown in Figs. 4-3, 4-4 and 4-5, respectively. 

The flow was initialized at freestream values everywhere except at the 
walls where a zero velocity condition was imposed. All features of the flow 
were then calculated by the GIM code. Convergence was assumed after 1500 
iterations when the sum of squares of the unsteady derivatives were decreasing 
by less than 1 percent. The turbulent viscosity was determined with the Baldwin- 
Lomax turbulence model and the laminar viscosity was held constant. The fr.ee- 

5 

stream Reynolds number per unit length was approximately 2 x 10 , 

These results will be used to form the initial conditions for the calcula- 
tion of the remainder of the turbulent flow field. 
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GUI CODE GRID 



Fig. 4-2 - Computational Grid for the Leading Edge of the Turbulent 
Spillage Problem 
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Fig. 4-3 - Spillage Problem - Leading Edge Adiabatic Walls 
(Turbulent) 
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Fig. 4-4 - Spillage Problem — Leading Edge Adiabatic Walls 
(Turbulent) 
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Fig. 4-5 - Spillage Problem — Leading Edge Adiabatic Walls 
(Turbulent) 
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5. CHEMISTRY MODEL DEVELOPMENT 


The original GIM algorithms were developed and the code written for 
an ideal gas only. The need for computing chemically reacting flow fields 
for application to engine combustion problems necessitates the development 
and implementation of a reacting gas capability in GIM. Most of the work 
during this contract has been limited to equilibrium considerations, but some 
attention has been given to finite rate models. The equilibrium work is dis- 
cussed first (Section 5.1) followed by a general formulation of a finite rate 
scheme for GIM. 

5.1 EQUILIBRIUM CHEMISTRY 

The basic premise of chemical equilibrium treatment in flow fields is 
that the local reaction rates are much faster than the flow residence time. 

The consequence is that the reactions proceed to completion at each point 
within the flow field. It has been shown in Ref 5-1 that for chemical equilibrium 
the thermochemistry computations can be uncoupled from the flowfield solu- 
tion. The thermochemistry data can then be communicated to the flowfield 
code via one of two methods. First, the flow state variables (P, p, T), 
thermodynamics (y, species distribution can be generated a priori 

for an isentropic expansion process from a given stagnation condition and 
the results stored in tabular form for use via a table look-up procedure. 

Second, the chemical equilibrium properties can be computed as an integral 
part of the flowfield solution by using the equilibrium calculation code as a 
subprogram. In either case the results are the same. Depending upon the 
size of the system, i.e., number of species and reactions involved, it has 
been found that, in general, the uncoupled approach saves computer time. 

A third approach is to simplify the reaction model to a complete reaction 
application. The approach is now discussed. 
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5.1»1 Complete Reaction Model 


Adding additional species equations to the GIM code and tracking each 
of the species present increases both the storage and computer time used in 
performing calculations for a reacting flow situation. As a first step, an 
effective global species conservation equation is solved for the mass fraction 
of the hydrogen fuel. This fuel mass fraction is the sum of the hydrogen com- 
ponent over all species; it is assumed that the remainder of the element mass 
fractions are those of nitrogen and oxygen (again summed over all species) in 
the proportions found in air. This approach allows the use of the two-gas species 
equation already in the GIM code. 

In the one and two ideal gas versions of the GIM code presently used, the 
thermodynamic and transport properties of the fluid have been computed locally 
whenever needed. This approach is efficient in the use of computer memory; 
however, when considering more complex chemistry models, this procedure 
can become too costly in the use of computer time. For this reason, the 
logic of the GIM code has been changed to perform the thermodynamic and 
transport property calculations only once per time step in a new fluid properties 
subroutine. For the sake of consistency and ease of use, the fluid property 
calculations for the one and two ideal gas cases have also been included in this 
subroutine, with the choice of paths through the subroutine being determined by 
an input indicator. 

The complete reaction model used to describe hydrogen-air chemistry 
assumes that all of the hydrogen fuel or oxygen present goes completely into 
water, depending upon whether there is an excess of oxygen or fuel, respectively. 
Following Ref, 5-2, reaction is not permitted when the volume fraction of 
hydrogen in the mixture is less than a threshold level of 0,04. Only the four 
species 02 > N 2 a^nd H 2 O are considered in the chemistry model des- 

cribed herein, and only when the amount of hydrogen present is less than the 
threshold level will both a*nd O 2 be present. The primary input to this 
phase of the complete reaction model is the mass fraction of hydrogen fuel 
(summed over all hydrogen-containing species) which is the dependent variable 
associated with the species equation. 
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After the species fractions have been determined, the temperature is 
found from the internal energy (obtained from the energy equation) using the 
enthalpy-temperature curve fits for the individual species from Ref. 5-3. 
Because these curve fits are higher order polynomials in temperature, an 
iterative procedure is required to do this. This iteration can be avoided as 
in Ref. 5-2 through the use of a quadratic equation to represent the enthalpy- 
temperature variation. This permits direct solution for temperature using 
the s tain dard inversion procedure for quadratic equations. 

Knowing the composition, temperature, and density (from the continuity 
equation), the mixture molecular weight can be obtained, and the pressure 
from the equation of state for a mixture of ideal gases. 

The transport properties are obtained by considering the fluid as an 
effective binary mixture for the purposes of computing a viscosity and thermal 
conductivity. Water, oxygen and nitrogen {H 2 O, O 2 and N 2 ) are lumped as one 
effective species and hydrogen (H 2 ) is the second species. Actual values of 
the transport properties are based on data from Refs. 5-4 and 5-5. Wilke's 
rule is used to obtain the mixture viscosity and Mason and Saxena's relation 
for the mixture conductivity. These relations were taken from Ref. 5-6. The 
transport coefficients needed for the species equation are obtained from con- 
stant Prandtl number and Lewis number. 

Since the speed of sound value used at various places in the GIM code 
need not be rigorously precise, an approximate value is computed from the 
perfect gas expression using a value for the ratio of specific heats of 1.4 
and a gas constant value obtained using the universal gas constant and the 
local mixture molecular weight. Frozen specific heat values are obtained 
using the curve fits in Ref. 5-2. 

The previous discussion has described how the complete reaction model 
composition and the associated thermodynamic and transport property data 
are computed. As noted earlier, the one and two ideal gas case fluid property 
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calculations have also been consolidated into the new fluid properties sub- 
routine. These ideal gas properties are computed as was done previously 
in the GIM code and are not described here. 

Initial checkout of the logic changes involved in implementing the fluid 
properties subroutine has been performed using a previously run one-ideal 
gas case, essentially reproducing the results obtained earlier. The hydrogen- 
air capability is being tested by running a flat-plate slot -injection test case 
(Ref. 5-7 as given in Ref. 5-2). Initially, results will be for non- reacting flow 
(obtained by setting the hydrogen fraction threshold level required for reaction 
to a value greater than one), and then for a reacting flow following the com- 
plete reaction model. 

5.1.2 Arbitrary Gas Full Equilibrium Model 

The mathematical formulation of a full chemical equilibrium calculation 
is well documented (Ref. 5-3) and a computer code is available to provide thermo- 
chemistry information. The computer code is fast computationally, reliable and 
accurate. Its output has been made compatible with our other flow field codes 
at Lockheed and are communicated to the codes via tabulated data. 

In the flow solution, the governing equations are cast in terms of inde- 
pendent variables such as velocity, entropy and total enthalpy. For isoener- 
getic isentropic flow there is no change in either the entropy or total enthalpy 
level during a flow field expansion process. Computationally, a combustion 
calculation (specified reactant formulation and pressure-temperature 
conditions) locates the appropriate isentrope where the flow field expansion is 
initiated (Fig. 5-1). Succeeding points on the isentrope are obtained by 
specifying one state variable (normally pressure) and the remaining properties 
computed from chemical equilibrium considerations and isentropic expansion 
calculations from reference conditions. 

Shock waves calculated in the discrete fashion are easily treated since 
the shock wave represents a discontinuity in the flow across which there is 
an increase in the entropy level as a result of loss in local total pressure, 
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This means simply that the expansion along an isentrope changes to a new isen- 
trope (Fig. 5-1) when the streamline of interest passes through a shock wave. 
To locate this isentrope, the change in entropy across the shock is computed 
from the shock relations. 

Chemical equilibrium is treated in the flowfield codes by generating 
tables of isentropic expansions a priori from specified reference conditions, 
where non- isentropic effects are obtained by expanding the isentropic relations. 
Most practical flow fields which are reacting have a non-constant energy 
field. This is treated in our equilibrium code as another parameter in a table 
look-up. An a priori knowledge of the range of energy distributions can be 
used to set up another entry in the table. This is commonly referred to as 
Oxidizer/Fuel (O/F) ratio (see Fig. 5-2). Utilization of the equilibrium 
routine then proceeds as described above. 

Shock capturing finite difference schemes are formulated by casting the 
set of governing partial equations in "conservative form" and then constructing 
the difference equation such that the conservative nature of the governing equa- 
tions is maintained. However, it is necessary to add artificial viscosity to the 
solution in the vicinity of shockwave for numerical stability. If too much arti- 
ficial viscosity is induced, the damping is too great and the jump conditions 
across the shock are computed incorrectly. Consequently, obtaining the 
correct thermodynamic properties through the shock depends on correctly 
constructing the shock (i.e., computing the correct "jump" conditions across 
the shock). In any event, the chemical equilibrium properties are obtained 
as discussed above, 

5.2 FINITE RATE MODEL DEVELOPMENT 

In flow problems where the gas may be considered in equilibrium 
(chemical and thermodynamic) at every point, two parameters are sufficient 
to define any of the other thermodynamic variables, either by assuming a 
perfect gas or by utilizing the results of chemical equilibrium calculations 
of the gases involved. If the gases cannot be considered to be in equilibrium. 


72 



Enthalpy, 


Composition 


® 1»®2 ^ 1*^2 


Input 


Note; OF is oxidizer -to-fuel ratio; 




j Chemical Equilibrium Code (CEC) | 
(Ref. 5-3) j 

I 1 


He,. Pc 




1 ^2 


Si 


Entropy, S 


Hc.t. P, 


H = total enthalpy 
P = pressure 


S = 
T = 
OF = 
A = 


entropy 
temperature 
oxidizer/fuel ratio 


Fig. 5-2 - Schematic of Equilibrium Thermochemistry Table Construction 
for Non-lsentropic and Non-Jsoenergetic Flow 




a point -by -point evaluation of its composition via integrating the individual 
species continuity equations is required. This subsection addresses a 
generalized procedure used to perform such calculations for an inviscid 
analysis This general procedure with addition of viscous terms will be 
used in the GIM code. 


A detailed description of the rate processes that occur in reacting 
flows requires that a myriad of mechanisms be considered to include all 
the possible chemical reactions of dissociation, formation, recombination. 
All of these, however, can be treated with a very general formalism. In 
the form usually quoted in chemical kinetics (Ref. 5-8) the phenomenological 
law of mass action states that the rate of a reaction is proportional to the 
product of the concentrations of the reactants. Thus, for a general reaction 
of the form 


N 

E^i 

i = l 


N 

A. 
1 1 


i = l 


^^.1) 


the net rate of production w^^ for any participating species for which the 
stoichiometric coefficients and p)' are not equal can then be written as 

N N 

w = k /7 (C ) -k n (C )^i ( 5 ^ 2 ) 

i=l ^ 1=1 ^ 

where is the species concentration defined as = p with p being the 
density and F^ being the species mol/mass ratio, respectively. Assuming 
small deviations from equilibrium, the forward and backward reaction rate 
constants, and kj^, respectively, can be related to the concentration equi- 
librium constant and to the pressure equilibrium constant as follows: 

N 

k, 2 (K - 

g- = = Kp ( 31T)‘ = ‘ (S.3) 
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The significance of the pressure equilibrium constant Kp is that it can be 
easily evaluated for any reaction using tabulated values of the equilibrium 
constant for formation from the elements. Values of are commonly tabu- 
lated in conjunction with specific heats, entropies and enthalpies as a function 
of temperature, and are available in general for most species. An equally 
convenient method exists for determining K from the change of free energy 

It 

accompanying the reaction, i.e,, 


K = exp(-AG/!KT) 
P 


(5.4). 


where AG is the change in free energy during the reaction process. Free 
energy values are also available for most species in tabular form. This is 
the method most commonly used to compute K . 


Using Eq. (5.3) to eliminate the backward rate constant from Eq. (5.2) 


the general production rate equation can be written as 


w 


N 

n (c.) 



N V\' 

77 (C.) 


(5.5), 


Finally, the net rate of production for any species participating in the 
reaction, either as a reactant or as a reaction product, is then given by 

= {v'/ - yl) w (5.6) 

Since most reaction systems involve a large number of simultaneous 
reactions, the net rate of production of species i usually equals a sum of 
terms. Thus, for an arbitrary number of M simultaneous reactions, all 
involving species i , 

M 

Wi = w k = 1 M (5.7) 
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where w. , 
1, k 

only. 


denotes the net rate of production of species i due to reaction K 


For reasons of computational speed and efficiency, a computer program 
can contain explicit expressions, as obtained from Eq. (5.5), for the most 
commonly encountered reaction mechanisms. Twelve types of reaction 
mechanisms are generally considered as possible contributors to the cal- 
culation of the net rate of production, w^ : 

Reaction type 


(1, 7) 

A + B 

7^ 

c 

+ D 

(2,8) 

A + B + M 


c 

+ M 

(3. 9) 

A + B 

— ^ 

c 

+ D + E 

(4, 10) 

A + B 

— 

c 


(5. 11) 

A + M 


c 

+ D + M 

(6, 12) 

A + M 

7^ 

c 

+ M 


Reaction types (7) through (12) correspond to reaction types (1) through 
(6), but proceed in the forward direction only. 

To reduce roundoff and truncation errors, the forward and backward 
rates for each reaction are computed separately. All contributions to the 
molar rate of production of a given species are then computed and added 
algebraically to form matrix coefficients. Since reaction types (7) through 
(12) proceed in the forward direction only, the second term on the right- 
hand side of Eq. (5.5) is disregarded in calculating the contributions to the 
coefficient matrix. 

In reactions (2), (5) and (6), as well as in (8), (11) and (12), M denotes 
a third body species which can be specified. For these reactions the situation 
often occurs where for various third bodies the respective rate constants 
differ only by a constant multiplier. These multipliers can be considered 
as third body efficiencies or weighting factors. If such a case is encountered, 


76 



the third body species mole mass ratio becomes effectively a fictitious 
mole mass ratio, consisting of the -weighted sum over all those species having 
a nonzero weighting factor, i.e.. 




'm = 4 -^ ^i ^M. 
1 1 


(5.9) 


where L are the weighting factors. 

The forward rate constant, k^, is generally a function of temperature. 

It is most often expressed in Arrhenius form. Again, for speed and efficiency 
in computation, the rate constants are divided into five types; 

Rate Constant Type 


(1) 

kf 

= A 

(2) 


= AT"^ 

(3) 

kf 

= A exp(B/ T) 

(4) 

kf 

= AT"’^ exp(B/ 5J T) 

(5) 

^f 

= AT“^ exp(B/5: T^) 


(5.10) 


The equations presented in this section provide a very general formalism 
for the evaluation of various rate processes. The specification of particular 
systems and associated rate constants is left to the program user to provide 
input data. 

Considering now the general species continuity equation 

pq. VC. = w (5.11) 

and making use of the foregoing discussion of the rate process we now proceed 
to describe a calculational technique for determining the individual species 
composition on a point -by -point basis. The description of this process is 
substantially simplified if Eq. (5,5) is specialized to a particular reaction 
type, say number (7) from Eq.(5.8) which is a one-way, two-body reaction. 
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A + B 


C + D 


(5.12) 


the net production rate for this process is 

w = - kj Fg (5.13) 

and the species continuity equation for species B then becomes 

Pq • P (5*14) 

which along streamlines becomes 

This equation can readily be solved using finite difference techniques 

employing explicit relationships, such as Euler or more sophisticated schemes, 

such as Predictor-Corrector. The step size for integrating this equation however 

is severely limited by stability criteria. It can be seen from Eq. (5.15) that the 

rate of change of a species along the streamline becomes increasingly larger 

as the flow speed is slowed, the density increased, or for fast reaction rates. 

In many flow problems, combinations of slow speeds, high densities and 

fast reaction rates (i.e., quasi -equilibrium) are quite common and integration 

- 8 

step sizes so small (i.e,, < 10 meters) are encountered that the solution 
becomes impractical in terms of computation time. 


For this reason, the technique described in Ref. 5-9» based on a 
linearization of the production rates, is a good choice. Writing Eq. (5.15) 
in finite difference form over a streamline step from station n to n+1. 


F 


n+1 



k^Sp 




( 5 , 16 ) 


And evaluating all the species concentrations at the downstream point results 
in a set of simultaneous nonlinear algebraic equations. In order to solve 
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these equations, we must then linearize the term ^Bn+l 

accomplished following Ref. 5-9, By expanding in terms of values at station 
n along with the increments over n to n+1 and neglecting products which are 

of second order, we can get the following expression using algebra. 


+^B ^A -^A ^'b 

n+l n+1 n n+1 n n+1 n 


Equation (5.l6) can now be written in its linearized form. Let C = AS k- p/q 



^B 

" ^B 

■ ^ K ^B ,, 

+ ^B ^A 


n+1 

n 

l n n+1 

n ] 

and 

""a ,, 

= ^A 

-^\^A , ^B 

+ F I 

B . , 


n + 1 

n 

1 n+1 

n+1 


n r 


(5.18) 


Equation (5.L8) can then be expressed in terms of a set of unknowns and 
calculable coefficients, C, Rewriting these we obtain 


where 


^B = V - (Fg -CFg (F^ ) 


n+1 


n + 1 


n 


n -^n+1 


n “n + 1 


= Qa - ^ ) - C F^ (F^ ) 


n 


B ' A ,, 
n n+1 


A '"B ,, 
n n+1 


Qj = Fj^ + C F, F, 
n n 


1 1 

n “'n 


(5.19) 


(5,20) 


^A ^ ^ C ^B ) + ^A ) ^B ,, = Qa 

n+:l n n n+1 n 


^A , ^B > ^ = "'a ' ^B ,, = ®B 

n+1 n n n + 1 n 


(5.21) 
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A matrix can now be formed using totally known information. 


l+CF^ 

n 

n 




«a" 

n 

_ n 

1 + C F . 

An 






The algebraic equations [A] [X] = [b] can then be solved for the unknown 

compositions F . , via an elimination technique. Although consuming 

n+1 n+1 

more time per integration step than an explicit formulation, the implicit tech- 
nique employed here is unconditionally stable permitting much larger step sizes. 

5.3 RECOMMENDATIONS 

Based on the development during the past year’s work, the following 
recommendations are offered. 


• The "complete reaction" equilibrium model should be 
checked out for hydrogen-air mixing in a duct. 

• The full, arbitrary gas, equilibrium model discussed in 
Section 5.1 should be completed and a solution compared 
to the complete reaction case. 

• Parallel to item 2 above, the full formal finite rate model 
should be coded in GIM. 

• A global finite rate system should be synthesized for 
hydrogen air to provide input for checkout of the finite 
rate coding. Comparison with the equilibrium cases 
should then be made. 
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6. INTERACTIVE GIM CODE INPUT PROGRAM ~ RUNGIM 


6.1 INTRODUCTION 

RUNGIM is an interactive program package designed to automate the 
preparation of GIM code program input and to effectively eliminate a 
percentage of the mistakes made by the novice as well as the experienced 
user. The prevention of gross errors will result in great savings of both 
engineering hours and computer time. In addition, the use of RUNGIM will 
reduce the required training period for a new user, make the GIM code easier 
to use, and significantly decrease the overall time required for setting up a 
given problem. 

RUNGIM combines the interactive features of FORTRAN Ebctended (FTN) 
Version 4,7 and Nos, 1.3 CYBER Controls Language (CCL) to efficiently auto- 
mate the GIM code input procedures. One of the most important of these 
features is the acceptance of list-directed input from a remote terminal. 

This type of input accepts integer, floating point, exponential, or alphanumeric 
data if separated by a comma or space. Thus the physical portion of the input 
phase becomes a function of user convenience and is quick, easy, and relatively 
free of possible error. To facilitate this type of input, RUNGIM is self-directing 
and prompts the user for input data. The "prompt" supplies the user with 
minimal necessary information about the input and then invites a response. 
RUNGIM contains appropriate default values for those inputs which have a 
frequently assumed value and checks user -supplied numerical inputs for 
proper format (i.e., integer or real) and range. If the user enters an improper 
input value or attempts a default where none is allowed, RUNGIM responds 
with appropriate messages and continues to prompt the user until an appropriate 
response is entered. Since many GIM code runs involve only minor changes 
in previously submitted input data, RUNGIM allows the user the access and 
XEDIT previously created data files. Finally, RUNGIM provides the user 
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with the ability to interrupt execution of the interactive program at any point 
by entering a special ** abort code.** The user is then able to resume execution 
at the point of interruption, restart at the beginning, or terminate execution 
completely. 

There are basically three types of data that must be generated for use 
in a GIM code run. These three types are job control cards (runstreams), 
UPDATE directives, and the actual program input data. Control card types 
are generally independent of the problem under consideration; therefore, 

RUNGIM can supply all necessary control card information except for account 
and charge numbers, core storage and run time requirements, and file name 
designations. By limiting user responsibility for knowledge of the job control 
language to such easily understood information most of the control card errors, 
especially those committed by novice users, are eliminated. 

Since modification of the GIM code is often necessary to accommodate 
special applications, RUNGIM allows the user to supply directives to the 
UPDATE processor from the remote terminal or from previously built files. 

This requires some user knowledge of both the UPDATE processor and the 
GIM code but provides a convenient means of assembling the UPDATE directives. 

The primary input to the GIM code is, of course, the program input data 
and this is quite naturally the greatest source of potential error for the user. 
RUNGIM prevents many of these errors by ensuring that the proper input loops 
are chosen and explained at each logical step according to the options requested. 
This eliminates the omission of data or improperly placed data. The use of 
the free formatted input scheme, which allows RUNGIM to properly format 
the information on a data file, eliminates the danger of badly formatted data 
stopping the job or producing erroneous results. 
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6.2 RUNGIM DETAILS AND FEATURES 


RUNGIM actually consists of seven separate programs which are exe- 
cuted interactively under the control of the main program, as illustrated in 
Fig. 6-1. These seven individual programs, GEOM, MATRIX, INTEG, GIMPLT, 
UPDATE, DYNDIM, and RNSTRM, perform the following tasks: 


GEOM 

MATRIX 

INTEG 

GIMPLT 

UPDATE 

DYNDIM 

RNSTRM 


supplies input data for the GIM code GEOM 
modules; 

supplies input data for the GIM code MATRIX 
modules (not yet operational); 

supplies input data for the GIM code INTEG 
modules 

supplies input data for the GIM code plotting 
program; (not yet operational) 

supplies modification directives to the CDC 
UPDATE processor; 

supplies input data for the GIM code dynamic 
dimensioning programs; and 

supplies the control cards required to execute 
any GIM code module on the NOS/CY203 computers. 


The basic layout of the individual programs is the same for each one 
and is illustrated for GEOM in Fig. 6-2. Each program has the capability 
of using previously created data files or interactively creating new one(s). 

The interactive FORTRAN portion of each program cein supply the control 
cards required to invoke the XEDIT'er and/or save the data files as permanent 
files. Note that the user can interrupt the execution of each program by enter- 
ing an "abort code" in response to any prompt. The user then has the option 
of resuming execution at the point of interruption, restarting the program, or 
terininating execution of that particular program. This feature could be ex- 
panded to allow branching to various points within a program during interrupt 
processing. 
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Program RUNGIM 


GEOM 

MATRIX 

IN TEG 

GIMPLT 

UPDATE 

DYNDIM 

RNSTRM 








Program GEOM 


Start 

GEOM 


or New^V, 
GEOM Data 


Create New 
GEOM Data File 
Interactively 


Access Old 
GEOM Data File 


“^dit GEOlvT 
^ Data 


Program 
Interrupt at 
Any Point 


Create Control 
Cards to Invoke 
XEDIT for GEOM 
Data File 


Continue 


Continu^ 

Restart 

Stop 


Save GEOM^ 
. Data File 


Create Control 
Cards to Save 
GEOM Data File 
as Perm. File 

■ — t ~ 

End FORTRAN X 
Program / 


Execute Control 
Card File to 
XEDIT and/or Save 
GEOM Data File 


Return to 
RUNGIM 
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The various features of RUNGIM are illustrated in the following examples 
of actual RUNGIM input and output. The program begins by displaying a banner 
message and prompts for display of instructions and execution code; 


tttttttttttttttttttttttttttxtttttttttttttttttttttttttttttttttttttttt 
BEGINNING RUNGIN VERSION DATE: 06/S6/8 

YOU WILL BE PROMPTED FOR INFORMATION REQUIRED TO CREATE OR ACCESS 
VARIOUS DATA FILES AND THE RUNSTREAM 

DO YOU UISH TO SEE INSTRUCTIONS BEFORE PROCEEDING (DEFAULT • NO) 

? 

SELECT RUNGIM MODULES BY THE FOLLOWING CODECS ): 

GEOM • 1 MATRIX • 3 INTEG - 3 GIMPLT - A 

UPDATE » 5 DVNDIM - 6 RNSTRM - 7 

ENTER MODULE CODECS ), ONE CODE/COLUMN 
? 1567 

YOU HAVE SELECTED THE FOLLOWING MODULES: 

GEOM UPDATE DVNDIM RNSTRM 

ENTER GO TO CONTINUE (DEFAULT • RESELECT MODULES) 

? GO 


If the user had responsed ^'^YES *^to the prompt for instructions, then RUNGIM 
would have displayed a description of the program and instructions for its use 
(see Section 3.2). In this particular case, the user has selected the GEOM, 
UPDATE, DYNDIM, and RNSTRM modules. RUNGIM begins execution with 
the GEOM module: 


tt NOW ENTERING GEOM MODULE tt 

YOU MAY EITHER ACCESS AN EXISTING (OLD) GEOMETRY DATA FILE, 
OR CREATE A NEW ONE 
ENTER OLD OR NEW (DEFAULT - NEW) 

? 

ENTER GEOMETRY DATA FILE NAME (DEFAULT ■ STRINP) 

? GEOMDAT 
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The user has selected creation of a new geometry input data file to be 
named "GEOMDAT." GEOM then begins to prompt for various input items: 


tt BEGINNING CREATION OF GEONETRV DATA FILE GEONOAT tt 
VOU UILL BE PROMPTED FOR INPUT 

ENTER HEADER — ANY ALPHANUNERIC INFORMATION MAY BE USED (DEFAULT - BLANK) 

? GEOMETRY DATA — FILE GEOMDAT 

SPECIFY NUMBER OF ZONES INTO UHICH THE FULL DOMAIN IS DIVIDED 
ENTER NZONES — NO LIMIT (DEFAULT • 1) 

7 a 

SPECIFY PROBLEM DIMENSIONALITY ~ ID, 3D, OR 3D 
ENTER 1, a, OR 3 (DEFAULT - 3) 

? 

SPECIFY ONE-STEP OR TUO-STEP TIME INTEGRATION SCHEME 
ENTER 1 OR a (DEFAULT - 3) 

•> 

SPECIFY GENERATION OF BOTH GRID AND FINITE DIFFERENCE MATRICES (0) 

OR GRID GENERATION ONLY (1) 

ENTER e OR 1 (DEFAULT - 0) 

? 

SPECIFY DO NOT MATE ZONES (0) OR MATE ZONES (1) 

ENTER e OR 1 (DEFAULT • 1) 

? 

SPECIFY ELLIPTIC RUN (0) OR QP RUN(l) 

ENTER 0 OR 1 (DEFAULT - 0) 

? 

SPECIFY NO DEBUG OUTPUT (0) OR PRINT DEBUG OUTPUT (1) 

ENTER 0 OR 1 (DEFAULT - 0) 

? 

SPECIFY NO PRINTOUT OF EACH ELEMENT MATRIX (0) 

OR PRINT EVERY NTH ELEMENT MATRIX (N) 

ENTER 0 OR N (DEFAULT - 0) 

7 

SPECIFY GRID POINT PRINTOUT FOR BOUNDARY NODES ONLY (0) 

OR FOR EVERY NTH NODE (N) 

ENTER 0 OR N (DEFAULT - 0) 


Note the extensive use of defaults for input items with commonly assumed 
values. This is further illustrated below: 
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SPECIFY FORWARD (F) OR BACKWARD (B) FINITE DIFFERENCES 


FOR EACH 
ENTER F 

TIME 
OR B 

STEP AND EACH COORDINATE DIRECTION 
FOR* 

STEP 

► 

1 

— 

X 

DIRECTION 

(DEFAULT 

- F) 

STEP 

1 

— 

Y 

DIRECTION 

(DEFAULT 

- F) 

STEP 

2 

— 

X 

DIRECTION 

(DEFAULT 

- B) 

STEP 

► 

2 

— 

Y 

DIRECTION 

(DEFAULT 

- B) 


BEGIN INPUT DATA FOR ZONE 1 

SPECIFY NUMBER OF SECTIONS WITHIN ZONE 1 
ENTER NSECTS — NO LIMIT (DEFAULT • 1) 

? 

SPECIFY NUMBER OF SEGMENTS IN EDGE 1 
ENTER NSEGS — 1 TO S (DEFAULT - 1) 

? 

SPECIFY EDGE SHAPE FUNCTION FOR EDGE 1 
USE THE FOLLOWING CODES* 


LINEAR 

• 1 

TRIG FUNCTION 


CIRCULAR ARC 

- 2 

OF X 

• 5 

CONIC SECTION 

- 3 

TRIG FUNCTION 


HELICAL ARC 

- 4 

OF THETA 

- 6 


SPECIAL FUNCTION • 7 


ENTER NSHAPE — 1 

TO 7 

(DEFAULT • 1) 



SPECIFY NUMBER OF SEGMENTS IN EDGE 3 
ENTER NSEGS — 1 TO 5 (DEFAULT - 1) 

? 

SPECIFY EDGE SHAPE FUNCTION FOR EDGE 3 
ENTER NSHAPE — 1 TO 7 (DEFAULT ■ 1) 

? 2 


SPECIFY NUMBER OF SEGMENTS IN EDGE 3 
ENTER NSEGS — 1 TO 5 (DEFAULT - 1) 

? 

SPECIFY EDGE SHAPE FUNCTION FOR EDGE 3 
ENTER NSHAPE — 1 TO 7 (DEFAULT - 1) 

? 

SPECIFY NUMBER OF SEGMENTS IN EDGE 4 
ENTER NSEGS — 1 TO 5 (DEFAULT • 1) 

? 

SPECIFY EDGE SHAPE FUNCTION FOR EDGE 4 
ENTER NSHAPE — 1 TO 7 (DEFAULT - 1) 

? 2 



Some input prompts display tables of appropriate values for the initial 
prompt only; subsequent prompts for the same input do not display the table. 
The edge shape function prompt shown above and the boundary condition indi- 
cator prompt shown below are examples of this ’’table on first prompt" pro- 
cedure. 


SPECIFY BOUNDARY CONDITION INDICATOR FOR EDGE 


USE THE FOLLOWING CODES t 
CONSTANT NODES ■ 0 
AXIS NODES - 1 
NO-SLIP/STAG- 
NATION NODES - 2 
CORNER NODES • 3 

INTERIOR NODES 
ENTER IBUL — 0 TO 9 (DEFAULT - 


FREE-SLIP/ 
TANGENCV 
SPECIAL CASE 
ONE-SIDED 
DIFFERENCES 
- 9 
9) 


4 

5-7 

8 


Subsequent prompts for NSHAPE and IBWL will not display the table of values 
in order to expedite the input process. 


The following sequence illustrates the program' s ability to respond to 
improper input: 


SPECIFY NUMBER OF NODES IN THE ETAl DIRECTION 
ENTER NNOD — 2 TO 100 (NO DEFAULT) 

? 

NO DEFAULT FOR NNOD 

SPECIFY NUMBER OF NODES IN THE ETAl DIRECTION 
ENTER NNOD — 2 TO 100 (NO DEFAULT) 

? 0 

INPUT ERROR — NNOD 

PLEASE OBSERUEt 2 .LE* NNOD .LE« 100 

SPECIFY NUMBER OF NODES IN THE ETAl DIRECTION 
ENTER NNOD — 2 TO 100 (NO DEFAULT) 

? 99.95 

99. < ERROR, RETYPE RECORD AT THIS FIELD 

? 99 
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At the first request for the number of nodes (NNOD), the user attempts to 
default. GEOM responds by indicating that there is no default for NNOD and 
repeats the prompt. This time the user enters an illegal value (0). GEOM 
indicates an error, displays proper bounds for NNOD, and repeats the prompt. 
The user next enters a value of improper format (real value instead of the 
integer value implied by NNOD). Again GEOM indicates the error and re- 
prompts. This time the user enters a legal value in the format consistent 
with the input variable and the program continues. 


The virtues of the list-directed input mode are illustrated below where 
several values are requested simultaneously; 


INPUT POINT COORDINATES, FLOU ANGLES, AND SEGMENT EXTREMALS (IF ANV 
OBSERUE THE FOLLOUING SEQUENCE FOR 
COORDINATES i 

PTl • X COORDINATE OF POINT 1 

PT2 ■ V COORDINATE OF POINT 1 

PT3 • Z COORDINATE OF POINT 1 

PT4 • FLOU ANGLE IN THE X-V PLANE AT POINT 1 

PT5 ■ FLOU ANGLE IN THE X-2 PLANE AT POINT 1 

ETC., FOR EACH POINT 
INPUT COORDINATES OF POINT 1 
ENTER PT ARRAY FOR POINT 1 (5 UALUES/LINE) 

? 0.5 0.0 0.0 0,0 0,0 

INPUT COORDINATES OF POINT 2 
ENTER PT ARRAY FOR POINT 2 (5 UALUES/LINE) 

? 1.0 0.0 0.0 0,0 0,0 


Here several values are input on the same line without regard to column loca- 
tion. Either blanks or commas may be used to separate values. Default, 
range, and format checking and error processing are invoked for each value. 

GEOM continues prompting for and accepting input data for the grid 
generation program, executing all of the loops required to create a complete 
input data file. The beginning of each major loop in the program is indicated 
so that the user knows precisely what portion of the geometry the subsequent 
input describes; 
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BEGIN INPUT DATA FOR ZONE 


2 


SPECIFV NUI1BER OF SECTIONS UITHIN ZONE 2 
ENTER NSECTS — NO LIMIT (DEFAULT • 1) 

? 

SPECIFV NUMBER OF SEGMENTS IN EDGE 1 
ENTER NSEGS — 1 TO 5 (DEFAULT ■ 1) 


The next example illustrates how GEOM can access an existing geom- 
etry input data file and then edit it via XEDIT. 


NOU ENTERING GEOM MODULE tt 

VOU MAV EITHER ACCESS AN EXISTING (OLD) GEOMETRY DATA FILE, 

OR CREATE A NEU ONE 

ENTER OLD OR NEU (DEFAULT - NEU) 

? OLD 

ENTER GEOMETRY DATA FILE NAME/USER NUMBER (DEFAULT • t/LOGIN NUMBER) 
? GEOMDAT 

DO YOU UISH TO XEDIT FILE GEOMDAT (DEFAULT • NO) 

? YES 

DO YOU UISH TO SAUE FILE GEOMDAT (DEFAULT • NO) 

? 


PREPARE TO XEDIT FILE GEOMDAT 
XEDIT 2.1.7 
?? Pt 


2-D SUPERSONIC SOURCE FLOU CASE — A TEST FOR RUNGIM 


2 

2 

0 

1 




1 0 

1 






F 


B 

B 




1 







2 

1 

2 





8 

4 

0 

0 

0 

1 


11 

e 

0 

0 

0 



o.e 


0.0 


0.0 

0.0 


e.e 


0.0 


0.0 

0.0 


.5 


0.0 


0.0 

0.0 

0.0 

1. 


0.0 


0.0 

0.0 

0*0 

96593 

» i 

25882 


0.0 

15. 

0.0 

48296 

♦ 

12941 


0.0 

15. 

0.0 


/EOR 

END OF FILE 
77 
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Here the data file can be modified using XEDIT to affect any changes desired. 
After editing, the file could be saved; however, in this instance the user relied 
on the default "NO SAVE" option so that the file will remain as a local file only. 

The interrupt capability of RUNGIM is illustrated in the next example. 
After several prompts, the user decides that he really wanted to use an existing 
data file instead of creating a new one. In response to the next prompt, the 
user enters the "abort code" (-999 in this version): 


tt NOU ENTERING GEON NODULE tt 

YOU NAY EITHER ACCESS AN EXISTING (OLD) GEONETRY DATA FILE, 

OR CREATE A NEU ONE 

ENTER OLD OR NEW (DEFAULT - NEU) 

? 

ENTER GEOMETRY DATA FILE NAME (DEFAULT • STRINP) 

? GEONDAT 

tt BEGINNING CREATION OF GEONETRY DATA FILE GEONDAT tt 
YOU WILL BE PROMPTED FOR INPUT 

ENTER HEADER -- ANY ALPHANUMERIC INFORMATION HAY BE USED (DEFAULT • BLANK 
? I HAVE MADE A MISTAKE — I REALLY WANT TO USE AN OLD FILE 
SPECIFY NUMBER OF ZONES INTO WHICH THE FULL DOMAIN IS DIVIDED 
ENTER NZONES — NO LIMIT (DEFAULT - 1) 

? -999 


tt EXECUTION INTERRUPTED ~ MODULE GEOM tt 
YOU HAY CHOOSE TOt (1) CONTINUE EXECUTION OF GEOM 

(2) RESTART EXECUTION OF GEOM 

(3) TERMINATE EXECUTION OF GEOM 


ENTER 1, 2, OR 3 
? 2 


tt NOU ENTERING GEOM MODULE tt 

YOU MAY EITHER ACCESS AN EXISTING (OLD) GEOMETRY DATA FILE, 
OR CREATE A NEW ONE 
ENTER OLD OR NEU (DEFAULT - NEU) 
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GEOM responds by interrupting execution and giving the user several options, 
i.e., continuing execution, restarting, or terminating the program. The user 
selects a restart and the program restarts at the beginning of GEOM. 

All of the preceding examples have illustrated the GEOM program. 
However, the others are similar in procedures and features though they per- 
form different tasks. Several examples are given on the following pages to 
demonstrate the features and capabilities of the various RUNGIM modules. 

The sample problem used in these examples is the two-dimensional super- 
sonic source flow case of Ref. 1 (Section 8.1, Ref. 1). 

6.3 RUNGIM EXAMPLES 
6.3.1 Access and Execution 

The RUNGIM program package is available to any remote terminal user 
and may be accessed and executed using the following commands (in batch mode): 


/GET ( RUNGIf1/UN-750978C ) 
/RUNGIM. 


6.3.2 RUNGIM Instructions 

The following are the instructions issued by RUNGIM if requested by the 

user: 
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/RUNGin 


tttttttttttttttttttttttttttttttttttxtttttttttttttttttttttttttttttttttt 
BEGINNING RUNGIN VERSION DATE: 06/26/81 

VOU UILL BE PROMPTED FOR INFORMATION REQUIRED TO CREATE OR ACCESS 
VARIOUS DATA FILES AND THE RUNSTREAH 

DO VOU UISH TO SEE INSTRUCTIONS BEFORE PROCEEDING (DEFAULT • NO) 

? YES 


tt RUNGIM INSTRUCTIONS ** 


PLEASE SELECT AN ITEM FROM THE FOLLOUING MENU: 

1. GENERAL DESCRIPTION OF RUNGIM 

2. CATALOG OF RUNGIM MODULES 

3. FREE-FORMAT INPUT DESCRIPTION 
A, ERROR PROCESSING (ABORT MODE) 

5. RUNGIM DEFAULTS 

6. COMPATIBILITY UITH CURRENT GIM CODE DECKS 

7. KNOUN ■BUGS' IN RUNGIM 

8. SUGGESTIONS AND COMMENTS 

9. ALL OF THE ABOVE 

10. RETURN TO RUNGIM 

ENTER ITEM NUMBER (DEFAULT - 10) 

? 9 
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GENERAL DESCRIPTION OF RUNGIN 

RUNGin IS A COLLECTION OF FORTRAN PROGRANS AND CCL PROCEDURES 
WHICH PERMIT THE USER TO INTERACTIUELY CREATE INPUT DATA FILES 
AND RUNSTREANS FOR THE UARIOUS GIN CODE NODULES. RUNGIN QUERIES 
THE USER FOR INPUT DATA UIA SELF-EXPLANATORY PROMPTS. AS THE USER 
ENTERS THE DATA, RUNGIN ARRANGES THE DATA IN THE PROPER SEQUENCE 
AND FORMAT AND WRITES IT OUT ONTO USER-NAMED FILES. ALL USER- 
SUPPLIED INPUT IS SCANNED FOR PROPER FORMAT (INTEGER, REAL, OR 
ALPHANUMERIC) AND RANGE; THE USER IS NOTIFIED OF ERRORS AND 
INUITED TO REENTER CORRECT DATA. AFTER CREATION OF A DATA FILE 
OR RUNSTREAM, THE USER HAS THE OPTION OF 'XEDIT'ING THE FILE 
AND/OR ■SAUE'ING IT. 

ENTER CARRIAGE RETURN TO CONTINUE 

? 

CATALOG OF RUNGIN MODULES 

RUNGIN HAS THE FOLLOWING MODULES: 

1. GEOM — CREATES INPUT DATA FILE FOR GEOMETRY MODULES 

2. MATRIX — CREATES INPUT DATA FILE FOR MATRIX MODULES (NA) 

3. INTEG — CREATES INPUT DATA FILE FOR INTEGRATOR MODULES 

4. GIMPLT — CREATES INPUT DATA FILE FOR PLOT MODULES (NA) 

5. UPDATE — CREATES DIRECTIUE FILE FOR THE UPDATE PROCESSOR 

6. DYNDIM — CREATES INPUT DATA FILE FOR THE DYNAMIC DIMENSIONERS 

7. RNSTRM — CREATES RUNSTREAMS FOR ALL GIM CODE MODULES 

ENTER CARRIAGE RETURN TO "CONTINUE 

? 
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FREE-FORWAT IHPUT DESCRIPTION 


RUNGIM USES FREE-FORMAT (LIST-DIRECTED) INPUT SO THAT THE USER 
NEED NOT BE CONCERNED ABOUT COLUMN NUMBERS OR SPACING. RUNGIM 
PROMPTS ARE OF THE SAME TYPE (INTEGER OR REAL) AS THE NUMERIC 
DATA TO BE ENTERRED. THUS THE USER MIGHT RESPOND AS FOLLOUSt 

. . . ENTER IDIM? 2 (INTEGER INPUT FOR INTEGER PROMPT) 

. . . ENTER EMU? 0.1 (REAL INPUT FOR REAL PROMPT) 

REAL NUMBERS MAY BE ENTERRED AS DECIMAL NUMBERS OR UITH E-TYPE 
NOTATION* E.G.* 0.1 OR l.E-1. ALPHANUMERIC INPUT DOES NOT 

REQUIRE ENCLOSING QUOTES OR HOLLERITH COUNT - ALL THIS IS HANDLED 
AUTOMATICALLY BY RUNGIM. 

ENTER CARRIAGE RETURN TO CONTINUE 

? 


ERROR PROCESSING (ABORT MODE) 

AT ANY TINE DURING THE USE OF RUNGIM* THE USER NAY ELECT TO HALT 
EXECUTION OF THE PROGRAM BYE ENTERRINQ THE ABORT CODEt -999 
RUNGIM UILL THEN SUSPEND EXECUTION AND QUERY THE USER AS TO HOU 
TO PROCEED. 

ENTER CARRIAGE RETURN TO CONTINUE 

? 
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RUNGin DEFAULTS 


MANV INPUT ITEMS HAUE BUILT-IN DEFAULT UALUES FOR COMMONLY SELECTED 
VALUES; THESE DEFAULTS MAY BE SELECTED BY TYPING A CARRIAGE RETURN 
IN RESPONSE TO A PROMPT. DEFAULT UALUES ARE INDICATED IN THE 
PROMPT BY THE MESSAGE* (DEFAULT ■ 0 . 1 ) (FOR EXAMPLE). 

SOME PROMPTS HAUE NO DEFAULT; THE USER MUST ENTER SOME VALUE. 

RUNGIM UILL NOT PERMIT THE USER TO CIRCUMVENT NO-DEFAULT ITEMS. 


ENTER CARRIAGE RETURN TO CONTINUE 

? 


COMPATIBILITY UITH CURRENT GIM CODE DECKS 

RUNGIM MODULES ARE COMATIBLE UITH CURRENT GIM CODE DECKS AS 
INDICATED BELOU* 

GEOM — GEOMC 

INTEG — INTEGD & INTEGE 

UPDATE — ALL VERSIONS 

DYNDIM — DYNMATC/DYNDIMC 

RNSTRM — VERSIONS C AND D OF ALL DECKS 


ENTER CARRIAGE RETURN TO CONTINUE 

? 


KNOUN "BUGS' IN RUNGIM 

THE FOLLOUING ARE KNOUN 'BUGS' IN THE RUNGIM MODULES* 

1. GEOM« UPDATE^ DYNDIM, AND RNSTRM HAVE NO KNOUN BUGS. 

2. INTEG IS A NEU MODULE AND MAY BE ERROR PRONE INITIALLY 

QP ROUTES MAY BE PARTICULARLY VULNERABLE. 

ENTER CARRIAGE RETURN TO CONTINUE 

? 
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SUGGESTIONS AND CONMENTS 


1. READ RUNGin PROMPTS CAREFULLY AND FOLLOW DIRECTIONS EXPLICITLY. 
DO NOT ANTICIPATE PROMPTSj UAIT FOR THE PROMPT AND THEN ENTER 
THE DATA REQUESTED. 

a. SOME PROMPTS REQUEST MORE THAN ONE PIECE OF DATA TO BE INPUT 
PER LINE. ENTER THE EXACT NUMBER OF DATA ITEMS REQUESTED, 

NO MORE, NO LESS. MULTIPLE INPUT ITEMS ON ONE LINE MAY BE 
SEPARATED BY BLANKS OR COMMAS. DO NOT USE SLASHES (/). 

3. REFER ALL QUESTIONS, COMMENTS, SUGGESTIONS, AND KNOWN BUGS TO: 

MICHAEL ROBINSON 
LOCKHEED - HUNTSUILLE 
305-837-1800 EXT-384 

ENTER CARRIAGE RETURN TO CONTINUE 

? 


DO YOU WISH TO REUIEU ANY ITEMS? 
ENTER YES OR NO (DEFAULT • NO) 

? 


SELECT RUNGIM MODULES BY THE FOLLOWING CODECS ): 

GEOM - 1 MATRIX ■ 2 INTEG ■ 3 GIMPLT - 4 

UPDATE • 5 DYNDIM • 6 RNSTRM • 7 

ENTER MODULE CODECS ), ONE CODE/COLUMN 

? 
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6.3.3 Supersonic Source Flow (2-D) — GEOM 

The following illustrates a complete RUNGIM setup of the input data 
and runstream required to execute the GIM code GEOM module for the super- 
sonic source flow example (Section 8.1, Ref. 6-1). 


/RUNGIM. 


tttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt 
BEGINNING RUNGIM UERSION DATEt 06/86/81 

VOU WILL BE PROMPTED FOR INFORMATION REQUIRED TO CREATE OR ACCESS 
VARIOUS DATA FILES AND THE RUNSTREAM 

DO YOU UISH TO SEE INSTRUCTIONS BEFORE PROCEEDING (DEFAULT - NO) 

? 

SELECT RUNGIM MODULES BY THE FOLLOUING CODECS)! 

GEOM • 1 MATRIX • 8 INTEG - 3 GIMPLT - 4 

UPDATE ■ 9 DYNDIM • 6 RNSTRM • 7 

ENTER MODULE CODECS). ONE CODE/COLUMN 
? 1967 

YOU HAVE SELECTED THE FOLLOUING MODULES! 

GEOM UPDATE DYNDIM RNSTRM 

ENTER GO TO CONTINUE (DEFAULT • RESELECT MODULES) 

? GO 


tt NOU ENTERING GEOM MODULE 

YOU MAY EITHER ACCESS AN EXISTING (OLD) GEOMETRY DATA FILE. 
OR CREATE A NEU ONE 
ENTER OLD OR NEU (DEFAULT - NEU) 

? 

ENTER GEOMETRY DATA FILE NAME (DEFAULT ■ STRINP) 

7 GEOMDAT 
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t* BEGINNING CREATION OF GEONETRV DATA FILE GEONDAT tt 
YOU UILL BE PROMPTED FOR INPUT 

ENTER HEADER -- ANY ALPHANUMERIC INFORMATION HAY BE USED (DEFAULT - BLANK) 

? e-D SUPERSONIC SOURCE FLOW CASE — A TEST FOR RUNGIM 

SPECIFY NUMBER OF ZONES INTO WHICH THE FULL DOMAIN IS DIVIDED 
ENTER NZONES — NO LIMIT (DEFAULT • 1) 

? 

SPECIFY PROBLEM DIMENSIONALITY — ID, SD, OR 3D 
ENTER 1, a, OR 3 (DEFAULT - 2) 

? 

SPECIFY ONE-STEP OR TUO-STEP TIME INTEGRATION SCHEME 
ENTER 1 OR a (DEFAULT • 3) 

? 

SPECIFY GENERATION OF BOTH GRID AND FINITE DIFFERENCE MATRICES (0) 

OR GRID GENERATION ONLY (1) 

ENTER 0 OR 1 (DEFAULT • 0) 

? 

SPECIFY DO NOT MATE ZONES (0) OR MATE ZONES (1) 

ENTER 0 OR 1 (DEFAULT • 1) 

? 

SPECIFY ELLIPTIC RUN (0) OR (JP RUN(l) 

ENTER 0 OR 1 (DEFAULT • 0) 

? 

SPECIFY NO DEBUG OUTPUT (0) OR PRINT DEBUG OUTPUT (1) 

ENTER 0 OR 1 (DEFAULT • 0) 

? 

SPECIFY NO PRINTOUT OF EACH ELEMENT MATRIX (0) 

OR PRINT EVERY NTH ELEMENT MATRIX (N) 

ENTER 0 OR N (DEFAULT • 0) 

? 

SPECIFY GRID POINT PRINTOUT FOR BOUNDARY NODES ONLY (0) 

OR FOR EVERY NTH NODE (N) 

ENTER 0 OR N (DEFAULT - 0) 

7 1 
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SPECIFY FORUARD (F) OR BACKUARD (B) FINITE DIFFERENCES 
FOR EACH TIME STEP AND EACH COORDINATE DIRECTION 
ENTER F OR B FOR* 

STEP 1 — X DIRECTION (DEFAULT ■ F) 

? 

STEP 1 ' V DIRECTION (DEFAULT - F) 

? 

STEP 2 — X DIRECTION (DEFAULT ■ B) 

9 

STEP 2 — V DIRECTION (DEFAULT - B) 


BEGIN INPUT DATA FOR ZONE 1 

SPECIFY NUNBER OF SECTIONS UITHIN ZONE 1 
ENTER NSECTS — NO LIMIT (DEFAULT • 1) 

7 

SPECIFY NUMBER OF SEGMENTS IN EDGE 1 
ENTER NSEGS — 1 TO 5 (DEFAULT - 1) 

? 

SPECIFY EDGE SHAPE FUNCTION FOR EDGE 1 
USE THE FOLLOWING CODES* 

LINEAR - 1 TRIG FUNCTION 

CIRCULAR ARC - 2 OF X -5 

CONIC SECTION - 3 TRIG FUNCTION 

HELICAL ARC -4 OF THETA ■ 6 

SPECIAL FUNCTION • 7 
ENTER NSHAPE — 1 TO 7 (DEFAULT - 1) 

? 

SPECIFY NUMBER OF SEGMENTS IN EDGE 2 
ENTER NSEGS — 1 TO 5 (DEFAULT • 1) 

? 

SPECIFY EDGE SHAPE FUNCTION FOR EDGE 2 
ENTER NSHAPE — 1 TO 7 (DEFAULT ■ 1) 

? 2 
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SPECIFY NUMBER OF SEGMENTS IN EDGE 3 
ENTER NSEGS — 1 TO 5 (DEFAULT • 1) 

? 

SPECIFY EDGE SHAPE FUNCTION FOR EDGE 3 
ENTER NSHAPE — 1 TO 7 (DEFAULT • 1) 

? 

SPECIFY NUMBER OF SEGMENTS IN EDGE 4 
ENTER NSEGS -- 1 TO 5 (DEFAULT » 1) 

? 

SPECIFY EDGE SHAPE FUNCTION FOR EDGE 4 
ENTER NSHAPE — 1 TO 7 (DEFAULT « 1) 

? a 

SPECIFY BOUNDARY CONDITION INDICATOR FOR EDGE 1 

USE THE FOLLOWING CODES: 

CONSTANT NODES - 0 FREE-SLIP/ 

AXIS NODES « 1 TANGENCY - 4 

NO-SLIP/STAG- SPECIAL CASE • 5-7 

NATION NODES • a ONE-SIDED 

CORNER NODES - 3 DIFFERENCES - 8 

INTERIOR NODES ■ 9 
ENTER IBUL — 0 TO 9 (DEFAULT ■ 9) 

? 4 

SPECIFY ZONE TO WHICH EDGE 1 IS MATED (0 FOR NO MATE) 
ENTER MATE — 0 TO 1 (DEFAULT • 0) 

? 

SPECIFY BOUNDARY CONDITION INDICATOR FOR EDGE 3 
ENTER IBUL — 0 TO 9 (DEFAULT ■ 9) 

? 8 

SPECIFY ZONE TO WHICH EDGE S IS MATED (0 FOR NO MATE) 
ENTER MATE — 0 TO 1 (DEFAULT - 0) 

? 

SPECIFY BOUNDARY CONDITION INDICATOR FOR EDGE 3 
ENTER IBWL — 0 TO 9 (DEFAULT • 9) 

? 4 

SPECIFY ZONE TO WHICH EDGE 3 IS MATED (0 FOR NO MATE) 
ENTER MATE — 0 TO 1 (DEFAULT • 0) 
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SPECIFY BOUNDARY CONDITION INDICATOR FOR EDGE 4 
ENTER IBUL — 0 TO 9 (DEFAULT • 9) 

? 0 

SPECIFY ZONE TO UHICH EDGE 4 IS HATED (0 FOR NO NATE) 

ENTER NATE — 0 TO 1 (DEFAULT - 0) 

? 

SPECIFY NODAL NUNBERING SEQUENCE ALONG COORDINATE DIRECTIONS 
USE THE FOLLOUING CODES t 

ETAS^ ETAl - 1 
ETAl, ETAS - S 

ENTER ITRAIN — 1 TO 3 (DEFAULT • 1) 

? 

SPECIFY NUMBER OF NODES IN THE ETAl DIRECTION 
ENTER NNOD — S TO 100 (NO DEFAULT) 

? 11 

SPECIFY NODAL DISTRIBUTION OPTION IN THE ETAl DIRECTION 
USE THE FOLLOUING CODES: 

UNIFORM SPACING • 0 

REDUCE ETAl SPACING - -1 

INCREASE ETAl SPACING • +1 

INPUT LOCATION OF NODES - NNOD 

ENTER ISTRCH — 0 TO 11 (DEFAULT - 0) 

? 

SPECIFY NUMBER OF NODES IN THE ETAS DIRECTION 
ENTER NNOD ~ 3 TO 100 (NO DEFAULT) 

? 11 

SPECIFY NODAL DISTRIBUTION OPTION IN THE ETAS DIRECTION 
ENTER ISTRCH — 0 TO 11 (DEFAULT ■ 0) 

? 

INPUT AC COEFFICIENTS FOR ZONE 1 
REFER TO TABLE 4-3 IN THE BLUE BOOK FOR CODES 
ENTER 4 AC UALUES FOR EDGE 8 
? 0.0 0.0 0.0 0.0 

ENTER 4 AC UALUES FOR EDGE 4 

7 0.0 0.0 0.0 0.0 
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INPUT POINT COORDINATES, FLOU ANGLES, AND SEGMENT EXTREMALS (IF ANY) 
OBSERUE THE FOLLOWING SEQUENCE FOR 
COORDINATES* 

PTl • X COORDINATE OF POINT 1 

PT2 - V COORDINATE OF POINT 1 

PT3 • Z COORDINATE OF POINT 1 

PT4 - FLOU ANGLE IN THE X-V PLANE AT POINT 1 

PT5 « FLOU ANGLE IN THE X-Z PLANE AT POINT 1 

ETC., FOR EACH POINT 
INPUT COORDINATES OF POINT 1 
ENTER PT ARRAY FOR POINT 1 (5 UALUES/LINE ) 

? 0.5 0.0 0.0 0.0 0.0 


INPUT COORDINATES OF POINT 2 
ENTER PT ARRAY FOR POINT 2 (5 UALUES/LINE) 


? 1.0 0.0 0.0 0.0 0.0 
INPUT COORDINATES OF POINT 3 
ENTER PT ARRAY FOR POINT 3 (5 UALUES/LINE) 

? 0.96593 0.25882 0.0 15.0 0.0 
INPUT COORDINATES OF POINT 4 
ENTER PT ARRAY FOR POINT 4 (5 UALUES/LINE) 

? 0.48296 0.12941 0.0 15.0 0.0 
DO YOU WISH TO XEDIT FILE GEOMDAT (DEFAULT 
? YES 


DO YOU WISH TO SAUE FILE GEOMDAT (DEFAULT - 


NO) 

NO) 


? 


GEOM COMPLETE — ESTIMATED NUMBER OF NODES - 121 

PREPARE TO XEDIT FILE GEOMDAT 
XEDIT 2.1.7 
?7 Pt 
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SUPERSONIC 

SOURCE 

FLOU 

CASE 

— A TEST FOR 

RUNGIM 

1 a 

a 0 

1 




e 0 

1 





F F 

B 

B 




1 1 






1 a 

1 a 





A 8 

4 0 

0 

0 

1 


LI 11 

0 0 

0 

0 



0.0 

0.0 


0.0 

0.0 


0.0 

0.0 


0.0 

0.0 


.5 

0.0 


0.0 

0.0 

0.0 

1. 

0.0 


0.0 

0.0 

0.0 

.96593 

.85888 


0.0 

15. 

0.0 

.48a96 

.18941 


0.0 

15. 

0.0 


/EOR 

“EOR— 

END OF FILE 
?? Q 

GEONDAT IS A LOCAL FILE 


tt NOU ENTERING UPDATE MODULE tt 

YOU MAY EITHER ACCESS AN EXISTING (OLD) UPDATE DIRECTIUE FILE, 
OR CREATE A NEU ONE 
ENTER OLD OR NEU (DEFAULT - NEU) 

7 

ENTER UPDATE DIRECTIUE FILE NAME (DEFAULT • UPDATE) 

? UPDSRC 
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tt BEGINNING CREATION OF UPDATE DIRECTIVE FILE UPDSRC tt 
A BLANK LINE UILL TERMINATE DIRECTIVE INPUT 
ENTER ONE DIRECTIVE PER LINE 
? XIDENT, SOURCE 
? t 
? t 

? * ANY UPDATE DIRECTIVES 

? * 

? * 

9 

DO YOU WISH TO XEDIT FILE UPDSRC (DEFAULT • NO) 

9 

DO YOU UISH TO SAVE FILE UPDSRC (DEFAULT - NO) 

? 


tt NOU ENTERING DYNMAT/DYNDIM MODULE ** 


SPECIFY TOTAL NUMBER OF NODES 
ENTER NN (DEFAULT - lEl ) 

? 

SPECIFY MAXIMUM NUMBER OF SPECIAL NODE TERMS 
ENTER NSP — 0 TO 131 (DEFAULT • 131) 

? 

SPECIFY MAXIMUM NUMBER OF NODES IN A CROSS-PLANE 
ENTER NXMAX (DEFAULT ■ 11) 

? 

SPECIFY MAXIMUM NUMBER OF NODES IN MATED PLANES 
ENTER NMATE (DEFAULT - 11) 

? 

SPECIFY NUMBER OF PLANES OF NODES SAVED PER RECORD 
ENTER NSAV (DEFAULT • 1) 

7 11 
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tt THE UPDATE DIRECTIUES AND DVNNAT INPUT DATA HAUE BEEN COMBINED tt 

THE COMBINED FILE NAME IS UPDSRC 

DO YOU UISH TO XEDIT FILE UPDSRC (DEFAULT • NO) 

? YES 

DO YOU UISH TO SAVE FILE UPDSRC (DEFAULT - NO) 

? 

PREPARE TO XEDIT FILE UPDSRC 
XEDIT 2.1.7 
?? Pt 

tIDENT, SOURCE 

t 

t 

t ANY UPDATE DIRECTIVES 

t 

t 

/EOR 

121 2 121 11 11 11 0 
/EOR 

END OF FILE 
?? Q 

UPDSRC IS A LOCAL FILE 


t* NOU ENTERING RNSTRM MODULE tt 

YOU MAY EITHER ACCESS AN EXISTING (OLD) RUNSTREAM, 
OR CREATE A NEU ONE 
ENTER OLD OR NEU (DEFAULT - NEU) 

? 

ENTER RUNSTREAM FILE NAME (DEFAULT - RNSTRM) 

7 RUNSRC 
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tt BEGINNING CREf^TION OF RUNSTRERN FILE RUNSRC tt 
VOU UILL BE PROnPTED FOR INPUT 


SPECIFY JOBNANE, FIELD LENGTH, AND TIME LIMIT 
ENTER JOBNAME/FL/TL (DEFAULT IS MODULE DEPENDENT) 

? RUNSRc/eeeoo/ieo 
ENTER USER NUMBER 
? 750978C 

ENTER CHARGE NUMBER 
? 101857 

ENTER GIM MODULE (S) TO BE RUN 
? GEOHC 

ENTER STAR FILE16 LENGTH (DEFAULT - 4) 

? 

ENTER STAR FILE18 LENGTH (DEFAULT ■ 12) 

? 

ENTER STAR FILE20 LENGTH (DEFAULT ■ 4) 

? 

ENTER STAR GEOMB CONTROLLEE FILE LENGTH 
? 600 

IF YOU UISH TO SAUE THE NOS DAVFILE, ENTER THE DAYFILE NAME 
? NOSDAY 

IF YOU UISH TO SAUE THE STAR DAYFILE, ENTER THE DAYFILE NAME 
? STRDAY 

DO YOU UISH TO XEDIT FILE RUNSRC (DEFAULT ■ NO) 

? YES 

DO YOU UISH TO SAUE FILE RUNSRC (DEFAULT ■ NO) 

? 

DO YOU UISH TO AUTOMATICALLY SUBMIT THIS JOB (DEFAULT • NO) 

? 

PREPARE TO XEDIT FILE RUNSRC 
XEDIT 2.1,7 
?7 Pt 


no 



/JOB 

/NOSEQ 

RUNSRc^cneeeee^Tiee. 

USER 

CHARGE (LRC ) 

GET(DYNnAT-DVNMATC/UN-838700C) 
6ET(OLDPL-GEOf1C/UN»838700C ) 
UPDATE(F,C-TAPE8) 

DVNHAT . 

COPVCF ( TAPE3, GEOtlS ) 

REUIND(GEOnS) 

RETURN ( OLDPL> TAPES, TAPES ) 

TOSTARC INPUT, GEOnS) 

DAVFILE(NOSDAV) 

REPLACE (NOSDAY) 

EXIT. 

DAYFILE( NOSDAY) 

REPLACE (NOSDAY) 

/EOR 

/READ,UPDSRC 

STORE 400SDS RUNSRC B 

STRSIDE,T100. 

FORTRAN ( I -GEOnS, B-GEONB, 0- LB ) 
REQUEST(FILE16/4,T-P ) 
REQUEST(FILE18/12,T-P) 
REQUEST(FILE19/ia,T-P) 
REQUEST(FILE20/4,T-P ) 
LOAD(GEOHB,CN-QEOnGO,600,GRLPALL- ) 
GEONGO. 

TOAS(Z-750978C,FILE18-BI,FILE20) 

DAYFILE(STRDAY) 

EXIT. 

TOAS(Z-750978C,FILE18-BI,FILE20) 

DAYFILE(STRDAY) 

EXIT. 

/EOR 

/READ,GEOnDAT 


/EOR 

REUERT. FROM EXRNST 
END OF FILE 
?? Q 

RUNSRC IS A LOCAL FILE 
REUERT. RUNGIM COMPLETE 
/ 
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6.3.4 Supersonic Source Flow (2-D) — INTEG 


The following illustrates a complete RUNGIM setup of the input data 
and runstream required to execute the GIM code INTEG module for the super- 
sonic source flow example (Section 8.1, Ref. 6-1). 

/RUNGIli. 


ttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttti 
BEGINNING RUNGIM UERSION DATE: 06/E6/81 

VOU UILL BE PROMPTED FOR INFORMATION REQUIRED TO CREATE OR ACCESS 
VARIOUS DATA FILES AND THE RUNSTREAM 

DO VOU UISH TO SEE INSTRUCTIONS BEFORE PROCEEDING (DEFAULT - NO) 

? 

SELECT RUNGIM MODULES BV THE FOLLOUING CODECS ): 

GEOM » 1 MATRIX ■ E INTEG - 3 GIMPLT ■ 4 

UPDATE - 5 DYNDIM - 6 RNSTRM ■ 7 

ENTER MODULE CODECS )> ONE CODE/COLUMN 
7 3567 

YOU HAVE SELECTED THE FOLLOUING MODULES: 

INTEG UPDATE DYNDIM RNSTRM 

ENTER GO TO CONTINUE (DEFAULT ■ RESELECT MODULES) 

? GO 


tt NOU ENTERING INTEG MODULE XX 

YOU MAY EITHER ACCESS AN EXISTING COLD) INTEGRATOR DATA FILE, 
OR CREATE A NEU ONE 
ENTER OLD OR NEU (DEFAULT - NEU) 

7 

ENTER INTEGRATOR DATA FILE NAME (DEFAULT • STRINP) 

7 INTDAT 
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tt BEGINNING CREATION OF INTEGRATOR DATA FILE INTDAT tt 
VOU UILL BE PROMPTED FOR INPUT 

ENTER ITITLE ANV ALPHANUMERIC INFORMATION MAV BE USED (DEFAULT - BLANK) 
? S-D SUPERSONIC SOURCE FLOW CASE — A TEST FOR RUNGIM 

SPECIFY PROBLEM DIMENSIONALITY — ID, 2D, OR 3D 
ENTER 1, 2, OR 3 (DEFAULT - 2) 

7 

SPECIFY TIME INTEGRATION SCHEME) 

(1) ONE-STEP ELLIPTIC 

(2) TUO-STEP ELLIPTIC 

(3) QUASI-PARABOLIC 
ENTER 1, 2, OR 3 (DEFAULT - 2) 

? 

SPECIFY MAXIMUM NUMBER OF TIME-STEP ITERATIONS 
ENTER ITMAX — NO LIMIT (NO DEFAULT) 

7 200 

SPECIFY FLOU FIELD PRINT CONTROL FLAG 

FLOU FIELD UILL BE PRINTED EVERY IPRNT-TH ITERATION 
ENTER IPRNT — 0 TO 200 (DEFAULT • 0) 

7 200 

SPECIFY FLOU FIELD SAVE CONTROL FLAG 

FLOU FIELD UILL BE SAVED (ON FILE 22) EVERY ITSAV-TH ITERATION 
ENTER ITSAV — 0 TO 200 (DEFAULT - 0) 

? 200 

SPECIFY COLD START RUN (0) OR RESTART RUN (N), 

UHERE N IS THE RECORD NUMBER OF THE STARTING DATA ON FILE 22 
ENTER ISTART (DEFAULT - 0) 
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SPECIFV PRINT OUTPUT TYPE* 

(1) ONE-LINE OUTPUT -- NODE, X, V, Z, RHO, U, U, U, Q, E, P, & IB 
(3) TUO-LINE OUTPUT — AS ABOVE PLUS GAM, CS, PHIX, PHIY, SOS, T, i N 
(3) SIMPLIFIED OUTPUT -- NODE, RHO, P, M, PHIX, PHIY, & IB 
ENTER lOTYPE — 1 TO 3 (DEFAULT - 1) 

? 3 

DO YOU UISH TO PRINT THE ZERO-TH ITERATION? 

ENTER YES OR NO (DEFAULT - YES) 

9 

SPECIFV ENGLISH UNITS ( 1 ) OR CQS NETRIC/DINENSIONLESS UNITS (3) 

ENTER lUNITS — 1 OR 3 (DEFAULT - 3) 

? 

SPECIFV STARTING ITERATION NUMBER 
ENTER ITSTRT (DEFAULT • 0) 

? 

SPECIFY TREATMENT OF VISCOSITY TERMS* 

(0) INVISCID EULER EQUATIONS 

(1) NUMERICAL VISCOSITY (DAMPING) ONLY 

(3) LAMINAR VISCOSITY AND DAMPING 

(3) BALDUIN/LOMAX TURBULENCE MODEL 
ENTER I Vise — 0 TO 3 (DEFAULT - 0) 

? 

SPECIFY INVISCID TYPE STARTING SOLUTIONS (0) OR VISCOUS/ 

BOUNDARY LAVER TYPE STARTING SOLUTION (1) FOR NOZZLE FLOUS 
ENTER IDIST — 0 OR 1 (DEFAULT - 0) 

? 

SPECIFY ONE IDEAL GAS (0) OR A TUO-GAS SYSTEM (1) 

ENTER ISPEC — 0 OR 1 (DEFAULT - 0) 

? 

SPECIFY SUPERSONIC-TVPE/ONE-SIDED DIFFERENCES (0) OR SUBSONIC-TYPE/ 

MASS BALANCING TECHNIQUE (1) FOR INFLOU AND OUTFLOU BOUNDARIES 
ENTER IDS — 0 OR 1 (DEFAULT • 0) 

? 
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SPECIFY PRINT SUMSQ INFORHATION EVERY N-TH ITERATION 
ENTER IPSQ (DEFAULT - 1) 

? 

SPECIFY SUtlSQ NORMALIZATION BY 1ST ITERATION VALUES (0), 

OR NO SUMSQ NORMALIZATION (1) 

ENTER INORM — 0 OR 1 (DEFAULT - 0) 

? 

SPECIFY CALCULATION OF SUMSQ FOR CONSERVED VARIABLES (0)« 

OR FOR PRIMITIVE VARIABLES (1) 

ENTER ISBSM — 0 OR 1 (DEFAULT -0) 

? 

SPECIFY TOTAL NUMBER OF NODES 
ENTER NN — NO LIMIT (NO DEFAULT) 

? 121 

SPECIFY NUMBER OF NODES IN THE X-COORDINATE DIRECTION (ZONE 1) 

ENTER NNX — 2 TO 100 (NO DEFAULT) 

? 11 

SPECIFY NODAL POINT INCREMENT IN THE X-COORDINATE DIRECTION (ZONE 1) 
ENTER NDX (NO DEFAULT) 

? 11 

SPECIFY NUMBER OF NODES IN THE Y-COORDINATE DIRECTION (ZONE 1) 

ENTER NNY — 2 TO 100 (NO DEFAULT) 

? 11 

SPECIFY NODAL POINT INCREMENT IN THE Y-COORDINATE DIRECTION (ZONE 1) 
ENTER NDY (NO DEFAULT) 

? 1 

SPECIFY NO SPECIAL TREATMENT OF SHARP EXPANSION CORNERS (0) 

OR NPM SHARP CORNERS ARE TO BE TREATED EXPLICITLY 
ENTER NPM — 0 TO 10 (DEFAULT - 0) 

? 

SPECIFY NUMBER OF ZONES (IF A MULTI-ZONE PROBLEM) 

ENTER KZONES — 1 TO 10 (DEFAULT - 1) 
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IF A CONSTANT TINE STEP IS DESIRED, 

ENTER DTINE (DEFAULT - VARIABLE TINE STEP SELECTED) 

? 0.5E-2 

SPECIFY CFL FACTOR — DTfllN - DTFAC t CFL 
ENTER DTFAC — 0.0 TO 1.0 (DEFAULT • 1.0) 

? 

SPECIFY TIME STEP UPDATE EVERY INCDT-TH ITERATION 
ENTER INCDT — 0 TO 200 (DEFAULT - 1) 

? 

SPECIFY EXPLICIT NETHOD (0) OR IMPLICIT METHOD (1) 

ENTER IMPL — 0 OR 1 (DEFAULT - 0) 

? 

SPECIFY PRANDTL NUMBER 
ENTER REALK (DEFAULT - 0.72) 

? 

SPECIFY SPECIFIC HEAT RATIO 
ENTER GAM (DEFAULT - 1.4) 

? 

SPECIFY MOLECULAR WEIGHT OF THE GAS 
ENTER UM (DEFAULT ■ GAS CONSTANT TO BE SPECIFIED NEXT) 

? 

SPECIFY IDEAL GAS CONSTANT FOR THIS GAS 
ENTER RK (NO DEFAULT) 

? 1.0 

SPECIFY FORWARD (F) OR BACKWARD (B) FINITE DIFFERENCES 

FOR EACH TIME STEP AND EACH COORDINATE DIRECTION 
ENTER F OR B FORJ 

STEP 1 — X DIRECTION (DEFAULT • F) 

? 

STEP 1 — Y DIRECTION (DEFAULT ■ F) 

? 

STEP 2 — X DIRECTION (DEFAULT - B) 

? 

STEP 2 — V DIRECTION (DEFAULT ■ B) 
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DO VOU UISH TO USE NOZZLE-FLOU TYPE INITIALIZATION? 
ENTER YES OR NO (DEFAULT • NO) 

? 


BEGIN FLOU FIELD INITIALIZATION — USE AS MANY SPECIFICATIONS AS REGUIRED 
SPECIFY FIRST NODE TO BE INITIALIZED BY SPECIFICATION 1 
ENTER NJ — 1 TO 121 (DEFAULT - TERMINATE FLOU FIELD INITIALIZATION) 

? 1 

SPECIFY NODE NUMBER INCREMENT 
ENTER INC (DEFAULT - 0) 

? 1 

SPECIFY TOTAL NUMBER OF NODES SET BY SPECIFICATION 1 
ENTER NTOT — 1 TO 121 (DEFAULT - 1) 

? 121 

SPECIFY INPUT OF VELOCITY COMPONENTS (0) OR TOTAL VELOCITY (1) 

ENTER ITAN — 0 OR 1 (DEFAULT ■ 1) 

? 

SPECIFY TYPE OF INITIALIZATION TO BE DONE BY SPECIFICATION 1 

(0) INPUT INITIAL CONDITIONS NEXT 

(1) USE NOZZLE FLOU INITIALIZATION OPTION (ELLIPTIC ONLY) 

(2) USE CONDITIONS ENTERRED FOR PREVIOUS SPECIFICATION 

(3) USE THE USERIP OPTION 
ENTER ITYPE ~ 0 TO 3 (DEFAULT • 0) 

? 

SPECIFY FLOU FIELD INITIAL CONDITIONS — SPECIFICATION 1 
ENTER THE FOLLOUINGt (NO DEFAULT) 

MASS DENSITY 
? 1.4 

TOTAL VELOCITY 

? 2.0 

STATIC PRESSURE 

? 1.0 


117 



00 


SPECIFY FIRST NODE TO BE INITIALIZED BY SPECIFICATION S 
ENTER NJ — 1 TO 131 (DEFAULT • TERMINATE FLOU FIELD INITIALI 

? 

BEGIN NODAL OUTPUT CONTROL SPECIFCATIONS — USE AS MANY AS REQUIR 
OUTPUT IS RESTRICTED TO 121 NODES 
SPECIFY FIRST NODE TO BE PRINTED BY SPECIFICATION 1 
ENTER N1 — 1 TO 121 (DEFAULT - TERMINATE FLOU FIELD OUTPUT S 
? 1 

SPECIFY NODE NUMBER INCREMENT 
ENTER IC (DEFAULT - 0) 

? 1 

SPECIFY TOTAL NUMBER OF NODES TO BE PRINTED BY SPECIFICATION 1 
ENTER NT — 1 TO 121 (DEFAULT • 1) 

? 121 

SPECIFY FIRST NODE TO BE PRINTED BY SPECIFICATION 2 
ENTER N1 — 1 TO 121 (DEFAULT - TERMINATE FLOU FIELD OUTPUT S 

? 

DO YOU UISH TO XEDIT FILE INTDAT (DEFAULT - NO) 

? YES 

DO YOU UISH TO SAUE FILE INTDAT (DEFAULT • NO) 

? 

PREPARE TO XEDIT FILE INTDAT 
XEDIT 2.1.7 
?? P* 
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END OF 

FILE 
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INTDAT 

IS 

A LOCAL FILE 









e e 


o.e 

0.0 


tt NOU ENTERING UPDATE NODULE tt 

YOU NAY EITHER ACCESS AN EXISTING (OLD) UPDATE DIRECTIUE FILE> 

OR CREATE A NEU ONE 

ENTER OLD OR NEU (DEFAULT - NEU) 

? OLD 

ENTER UPDATE DIRECTIUE FILE NANE/USER NUNBER (DEFAULT ■ */LOGIN NUMBER) 
? UPDSRC 

DO YOU UISH TO XEDIT FILE UPDSRC (DEFAULT • NO) 

? 



tt NOU ENTERING DVNHAT/DYNDIH MODULE tt 


ARE THE UPDATE DIRECTIVES AND DYNDIM INPUT DATA BOTH ON FILE UPDSRC 
ENTER YES OR NO (NO DEFAULT) 

? NO 

SPECIFY ELLIPTIC RUN (0) OR QP RUN (1) 

ENTER 0 OR 1 (DEFAULT « 0) 

? 

SPECIFY TOTAL NUMBER OF NODES 
ENTER NN (DEFAULT - 131) 

? 

SPECIFY MAXIMUM NUMBER OF SPECIAL NODE TERMS 
ENTER NSP “ 0 TO 121 (DEFAULT ■ 121) 

SPECIFY MAXIMUM NUMBER OF BOUNDARY NODES 
ENTER MNB — 0 TO 121 (DEFAULT ■ 121) 

? 

SPECIFY NUMBER OF RECORDS OF DATA ON FILE20 
ENTER IREC (DEFAULT • 1) 


tt THE UPDATE DIRECTIVES AND DYNDIM INPUT DATA HAVE BEEN COMBINED tt 

THE COMBINED FILE NAME IS UPDSRC 

DO YOU UISH TO XEDIT FILE UPDSRC (DEFAULT ■ NO) 

? YES 

DO YOU UISH TO SAVE FILE UPDSRC (DEFAULT • NO) 

? 

PREPARE TO XEDIT FILE UPDSRC 
XEDIT 2.1.7 
7? 


120 


XIDENT, SOURCE 

t 

t 

t ANV UPDATE DIRECTIVES 

t 

t 

/EOR 

lai 2 e 121 121 0 1 

/EOR 

END OF FILE 
?? Q 

UPDSRC IS A LOCAL FILE 


tt NOU ENTERING RNSTRfI MODULE ** 

VOU flAV EITHER ACCESS AN EXISTING (OLD) RUNSTREAM, 

OR CREATE A NEU ONE 

ENTER OLD OR NEU (DEFAULT - NEU) 

? 

ENTER RUNSTREAM FILE NAME (DEFAULT • RNSTRM) 

? RUNSRC 

tt BEGINNING CREATION OF RUNSTREAM FILE RUNSRC tt 
VOU UILL BE PROMPTED FOR INPUT 


SPECIFY JOBNAME, FIELD LENGTH, AND TIME LIMIT 
ENTER JOBNAME/FL/TL (DEFAULT IS MODULE DEPENDENT) 
? RUNSRC/60000/10e 
ENTER USER NUMBER 
?- xxxxxxc 
ENTER CHARGE NUMBER 
? xxxxxx 
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ENTER Gin nODULE(S) TO BE RUN 
? INTEGD 

ENTER NOS FILE17 NANE 
? TAPE17 

ENTER NOS FILEBO NANE 
? TAPE20 

ENTER STAR FILE2S LENGTH (DEFAULT - 3) 

? 

ENTER STAR INTGB BINARY FILE LENGTH (DEFAULT - 100) 

? 

ENTER STAR INTGB CONTROLLEE FILE LENGTH 
? 1000 

IF VOU UISH TO SAUE THE NOS DAYFILE, ENTER THE DAYFILE NAME 
? NOSDAY 

IF YOU UISH TO SAUE THE STAR DAYFILE, ENTER THE DAYFILE NAHE 
? STRDAY 

DO YOU UISH TO XEDIT FILE RUNSRC (DEFAULT - NO) 

? YES 

DO YOU UISH TO SAUE FILE RUNSRC (DEFAULT ■ NO) 

? 

DO YOU UISH TO AUTOHATICALLV SUBMIT THIS JOB (DEFAULT - NO) 

? 

PREPARE TO XEDIT FILE RUNSRC 
XEDIT 2.1.7 
?? Pt 
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/JOB 
/NOSEQ 

RUNSRC^cneeeeo^Tiee* 

USER 

CHARGE ( LRC ) 

GET(DVNDH1-DVNDIHC/UN-83870eC ) 

GET ( OLDPL- INTEGD/UN-838700C ) 

UPDATE(F,C-TAPE8) 

DVNDin. 

COPYCF ( TAPE3 , INTGS ) 

REUIND(INTGS) 

RETURN(0LDPL,TAPE3,TAPE8) 

ATTACH(FILE17-TAPE17) 

ATTACH(FILEa0-TAPE20) 

TOSTAR( INPUT, INTGS, FILE20,FILE17»BI//U ) 

DAVFILE(NOSDAV) 

REPLACE (NOSDAV) 

EXIT, 

DAVFILE( NOSDAV) 

REPLACE (NOSDAV) 

/EOR 

/READ,UPDSRC 

STORE 400SOS RUNSRC B 

STRSIDE,T100, 

REQUEST(FILE22/3,T-P) 

FORTRAN(I-INTGS,B-INTGB/100,O-LB) 

LOAD(INTGB,CN-INTEGO,1000 

,GRLP-»TMUEC,*PRIf1,»EBUF,»TURB,»UPROP,»UBUF,*BOUND,»STEP 

,GRLP-XXBUF1,QRLP-XXBUF2,QRLP-«TAU,XQPN0D,XSEC0RD 

,GRSP-XCNTRL,XUCOEF,l;TDATA,«UECP,XSQ,«Pn,XSUBSBC,XUSER 

, »STEPI , »QPCOH, »QPPRNT, »CVGCOri, »AXSVH 

,GR0L-«Q3nAP) 

INTEGO, 

T0AS(Z-FILE22) 

DAVFILE(STRDAV) 

EXIT. 

TOASCZ- FILE22) 


DAVFILE(STRDAV) 

EXIT. 

/EOR 

/READ,STRINP 

/EOR 

REVERT. FROM EXRNST 
END OF FILE 
?? Q 

RUNSRC IS A LOCAL FILE 
REVERT. RUNGin CONPLETE 
/ 
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6.4 CURRENT STATUS OF RUNGIM 


RUNGIM is operational and available to all remote terminal GIM code 
users who wish to use it. The MATRIX and GIMPLT modules of RUNGIM are 
not yet operational but will be in the near future. All RUNGIM users should 
occasionally check the "KNOWN BUGS" section of the instructions for any 
recently discovered problems in the code. Any problems, questions, comments, 
or suggestions should be directed to Michael Robinson, Lockheed-Huntsville, 
205-837-1800, ext 384. 
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7. CALCULATION OF THE FLOW FIELD ABOUT A WING-BODY 

CONFIGURATION 

7.1 INTRODUCTION 

Figure 7-1 depicts the wing body configuration for which the flow field 
was computed using the hyperbolic GIM/STAR computer code. The com- 
plete configuration is 0.355 m in length with a 70 deg swept wing. The upper 
wing surface is a flat plate aligned with the longitudinal center line of the body. 
The lower wing surface makes a 3-deg angle with the longitudinal center and 
a 8.712 deg normal to the wing leading edge. The forebody, which is 0.173 m 
in length, is described by a power law function. The cylindrical afterbody 
is 0.1822 m in length. The freestream conditions, which are shown in Fig. 

7-1, represent a flight altitude of approximately 30,480 m. Since these cal- 
culations represent the initial attempt to compute a wing -body configuration 
with GIM all computations were done for zero angle of attack. 


Section 7.2 presents a discussion of the development of the solution. An 
analysis of the final solution is given in Section 7.3. This discussion is divided 
into two parts: (1) the flow field over the forebody, and (2) the flow field over 

the afterbody. 

7.2 DEVELOPMENT OF THE SOLUTION 
7.2.1 Forebody Flow Field 

As a matter of convenience, the computational domain was constructed 
with the input boundary beginning .003 m aft of the forebody nose. The stag- 
nation region and leading edges could be treated with code, but it is not neces- 
sary here. We can easily use tables for these inputs. The outer boundary of 
the freestream was made to conform to the anticipated shape of the bow shock. 
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Since all computations were made for a pitch plane orientation (no yaw), a 
symmetry plane was employed. Figure 7-2 depicts the computational mesh 
on the body and symmetry plane. Figure 7-3 shows the body and outer sur- 
face grids while Fig, 7-4 presents the body, symmetry plane, and outer surface 
mesh points. The computational grid for a typical cross- plane is shown in 
Fig. 7-5. Equal nodal spacing was used in the cross-plane directions, 
while a mild stretching function was employed in the main flow direction in 
order to obtain the desired resolution of the bow shock in the nose region, 

7.2.2 Afterbody Flow Field 

The computational domain for the afterbody was constructed to adequately 
describe the wing without becoming prohibitatively large. As with the fore- 
body a symmetry plane was employed, and since base flow calculations are 
not feasible with a hyperbolic or parabolic method, the computational domain 
extends only to the wing trailing edge. A sharp wing leading edge was employed 
and the wing-body juncture has been modeled as a distinct three-dimensional 
external corner. Figures 7-6 and 7-7 present perspective views of the com- 
putational mesh on the wing and body surfaces. Figure 7-8 and 7-9 show the 
computational mesh on the upper and lower wing surfaces. Figure 7-10 is a 
side view of the mesh on the afterbody. The thickness of the wing can be 
seen in this view. In Fig. 7-11 a three-dimensional view of the computational 
mesh on the wing, the body, and the outer surfaces of the computational domain 
is shown. Figure 7-12 is a plot of the complete mesh on the final three com- 
putational planes. 

In order to accurately describe the wing, the upper and lower wing sur- 
face mesh points had to be at the same spatial location on the initial plane. 

This results in several dual points occuring at the intersection of the wing 
with the body (See Fig. 7-13). As Figs. 7-14 through 7-19 depict, the wing 
mesh spacing becomes larger as the trailing edge is approached. 

Since the computational mesh was very dense on the forward 20% of the 
wing, the finite difference analogs had to be exarnined closely in this area* 
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Because of the presence of the dual points, the analogs on the initial plane 
are somewhat suspect. However, since the flowfield values on the initial 
plane are input and not computed, thes'e analogs were not used directly in 
the calculations. The analogs on the second plane, which are influenced by 
the initial plane geometry, were examined carefully and are correct. 

7.3 RESULTS AND DISCUSSION 

7.3.1 Forebody Flow Field 

A 10830 -node computational grid with 361 nodes in a cross -plane was 
used to compute the flow field about the forebody. The solid boundaries were 
treated in an inviscid manner. The bow shock was attached to the nose and 
held fixed for the initial plane. For subsequent planes the shock location and 
jump conditions were determined from the governing equations. These com- 
putations required 210.0 CPU seconds on the CYBER 203 to obtain the final 
solution. Figure 7-20 shows the velocity vectors on the symmetry plane, 
and Fig. 7-21 provides a magnified view of the nose region. The bow shock 
lies close to the body near the nose, and dissipates as the flow expands back 
to freestream near the base. Figure 7-22 depicts the static pressure contours 
over the symmetry plane. Figure 7-23 shows the cross-plane static pressure 
contours on a plane that is .022 m from the nose. This figure depicts the 
near proximity of the bow shock to the missile body. Figure 7-24 presents 
the static pressures on the final cross-plane. The bow shock has dissipated 
by this time, and the flow has expanded back to freestream. 

7.3.2 Afterbody Flow Field 

The computational grid for the afterbody contained 9796 nodes with 3l6 
nodes in each cross plane. The solid boundaries with the exception of the 
wing leading edge and the wing-body juncture point were treated in an inviscid 
manner. The wing leading edges were initialized to inviscid post shock values 
and held fixed during the calculations. This was done for convenience only 
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since the downstream shock conditions can be determined from shock tables. 
Tangency conditions could also have been used; however, this would require 
additional mesh points to permit the formation of the shock. The wing- body 
juncture presented a unique problem and was dealt with in a special manner. 
Inviscid corner flow conditions were used at the wing-body juncture. As 
subsequent planes were computed the mesh spacing in the vicinity of the 
corner became larger and larger. The divergence of the mesh is due to the 
incompatibility of the grid structure over different components of aircraft 
configurations. The solution to this problem is to use grid imbedding tech- 
niques to achieve grid resolution where it is required; however, time and 
budget restrictions did not permit this refinement. Inviscid boundary con- 
ditions could normally be applied if the mesh were tight enough. Since this 
refinement was not possible, the corner flow field conditions were held fixed 
in time at the final values of the first five corner points. The first five points 
along the wing/body juncture appeared to be fine because they were not affected 
by the mesh divergence problem. The implementation of the above two "speciar* 
boundary conditions did not create any perturbations in the computed flow 
field. These computations required 334 CPU seconds on the CYBER 203 to 
obtain the final solution. For more details of the computational procedure 
utilized in the wing/body problem please refer to Ref. 7-1. 

Cross-plane velocity vector plots for planes 5, 10, 15, 20, 25 and the 
final plane are presented in Figs. 7-25 to 7-30, respectively. The fact that 
the velocity vectors on the lower wing surface seem to be coming directly 
out of the wing is an illusion created by the fact that only cross-plane com- 
ponents are plotted. Each of these vectors has a much larger component 
directed into the plane of the plots. These plots depict the presence of the 
formation of vortices particularly from the lower wing surface. 

Static pressure contours are presented for planes 5, 10, 15, 20, 25 and 
the final plane in Figs. 7-31 through 7-36, respectively. A high pressure 
region is produced on the lower wing surface while the upper surface remains 
virtually at freestream conditions. 
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Figures 7-37 and 7-38 show the velocity vectors and static pressure 
contours, respectively on the wing upper surface for the grid of Fig. 7-8. 
Figures 7-39 and 7-40, respectively, depict the same type of plots on the 
wing lower surface for the grid of Fig. 7-9. The velocity vectors and pres- 
sure contours are shown in Figs. 7-41 and 7-42, respectively, for the grid 
of Fig. 7-10. 

7.3.3 Comparison with Experimental Data 

Experimental data for this configuration was obtained from the NASA- 
Langley Research Center. 

Figures 7-43 and 7-44 present a comparison between GIM results and 
this experimental data for the wing/body configuration. As Fig. 7-43 depicts, 
the GIM results compare extremely well with the experimental data with GIM 
predicting slightly more over-expansion than the data. The spanwise pres- 
sure distribution in Fig. 7-44 shows good comparison between GIM and exper- 
iment, The slight deviation between the GIM predictions and the experimental 
data near the wing leading edge is due to the imposition of the fixed post shock 
boundary condition on the wing leading edge. 


SECTION 7 REFERENCE 


7.1. Xiques, K. E. et al., ** Computation of Three-Dimensional Inviscid 
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131 





Fig. 7-3 - For^ebody Computational Mesh (Body Surface and Outer 
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j'jg^7_4 - Forebody Computational Mesh (Body Surface, Symmetry 
Plane, and Inter Surface) 
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Fig. 7-5 - Computational Mesh for a Typical Cross-Plane 
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Fig. 7-6 - Computational Mesh on Wing-Body Surface 



Fig. 7-7 - Computational Mesh on Wing-Body Surface (aft view) 
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Fig. 7-8 - Computational Mesh on the Upper Wing Surface Viewed from Above 
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Fig. 7-9 


- Computational Meah on the Lower Wing Surface Viewed from Below 
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Pig. 7-10 - Computational Mesh on Afterbody Surface 
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Fig. 7-11 - Computational Mesh on the Wing-Body and Outer 
Surfaces (Aft View) 
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7-12 - Three-Dimensional Grid for Final Three Computation 
Planes (Aft View) 
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Fig. 7-13 - Afterbody Cross-Plane Grid, Initial Plane, X = 0.1730 m 
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Fig. 



^.14 - Afterbody Cross 


-Plane Grid, Plane 5, X 


= 0. 1973 m 
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Fig. 7-15 - Afterbody Cross-Plane Grid, Plane 10, X = 0,2276 


m 
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Fig. 7-16 - Afterbody Cross-Plane Grid, Plane 15, X - 0.2641 m 
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Fig. 7-17 - Afterbody Cross-Plane Grid, Plane 20, X 


0,2884 m 
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Fig. 7-18 - Afterbody Cross-Plane Grid, Plane 25, X = 0.3188 m 
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Fig. 



-19 - Afterbody Cross-Plane Grid, Final Plane, X - 0.3551 m 
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Fig. 7-20 - Forebody Velocity Vectors (Symmetry Plane) 
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Fig. 7-22 - Static Pressure Contours (Symmetry Plane) 
(Pq =1.0 lbf/ft2 = 0.04788 k Pa) 


153 



ID 


1 

S 


p/p 



0000 

0000 

0000 



Fig. 7-24 
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-Plane Static Pressure Contours (Final Plane) 
1,0 Ibf/ft^ = 0.04788 ,k Pa) 
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Fig. 7-35 - Static Pressure Contours for Plane 25, X = 0.3188 
(P^ = 1.0 Ibf/ft^ = 0.04788 k Pa) 
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Fig. 7-40 - Lower Wing Surface Static Pressure Contours 
(Pq = 1.0 lbf/ft2 = 0.04778 k Pa) 
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Fig. 7-43 - Surface Pressure Profile on Forebody 
(P = 24.3 lbf/ft2 = 1.1634 k Pa) 

(Zq = 1.0 ft = 0.3048 m) 



Spanwise Distance Along Wing 


Fig. 7-44 - Spanwise pressure Distribution 
Z = 0.2804 m 

(P = 24.3 Ibf/fr = 1.1634 k Pa) 
(Y°° = 1.0 ft = 0.3048 m) 
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